#sequencing data used to generat OTU tables via Qiime2 version 2018.11, can be found under NCBI SRA accession: PRJNA615043
getwd()
## [1] "/Users/martinylab/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020"
setwd("/Users/martinylab/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020")
library(BiodiversityR)
## Loading required package: tcltk
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
## Loading required package: vegan3d
## Loading required package: rgl
## Registered S3 methods overwritten by 'lme4':
##   method                          from
##   cooks.distance.influence.merMod car 
##   influence.merMod                car 
##   dfbeta.influence.merMod         car 
##   dfbetas.influence.merMod        car
## BiodiversityR 2.11-3: Use command BiodiversityRGUI() to launch the Graphical User Interface; 
## to see changes use BiodiversityRGUI(changeLog=TRUE, backward.compatibility.messages=TRUE)
library(biomformat)
library(car)
## Loading required package: carData
library(data.table)
library(doBy)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
## 
##     between, first, last
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(devtools)
## Loading required package: usethis
## 
## Attaching package: 'devtools'
## The following object is masked from 'package:permute':
## 
##     check
library(EcolUtils)
library(emmeans)
## 
## Attaching package: 'emmeans'
## The following object is masked from 'package:devtools':
## 
##     test
library(forcats)
library(gapminder)
library(ggplot2)
library(ggpubr)
library(lmerTest)
## Loading required package: lme4
## Loading required package: Matrix
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(lsmeans)
## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:dplyr':
## 
##     intersect, setdiff, union
## The following objects are masked from 'package:data.table':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(multcompView)
library(purrr)
## 
## Attaching package: 'purrr'
## The following object is masked from 'package:data.table':
## 
##     transpose
## The following object is masked from 'package:car':
## 
##     some
library(qiime2R)
library(readxl)
library(rcompanion)
library(Rmisc)
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following object is masked from 'package:purrr':
## 
##     compact
## The following object is masked from 'package:ggpubr':
## 
##     mutate
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
library(stringr)
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following objects are masked from 'package:Matrix':
## 
##     expand, pack, unpack
library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ tibble 3.0.1     ✓ readr  1.3.1
## ── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x plyr::arrange()          masks dplyr::arrange()
## x lubridate::as.difftime() masks base::as.difftime()
## x dplyr::between()         masks data.table::between()
## x readr::col_factor()      masks scales::col_factor()
## x plyr::compact()          masks purrr::compact()
## x plyr::count()            masks dplyr::count()
## x lubridate::date()        masks base::date()
## x scales::discard()        masks purrr::discard()
## x tidyr::expand()          masks Matrix::expand()
## x plyr::failwith()         masks dplyr::failwith()
## x dplyr::filter()          masks stats::filter()
## x dplyr::first()           masks data.table::first()
## x lubridate::hour()        masks data.table::hour()
## x plyr::id()               masks dplyr::id()
## x lubridate::intersect()   masks base::intersect()
## x lubridate::isoweek()     masks data.table::isoweek()
## x dplyr::lag()             masks stats::lag()
## x dplyr::last()            masks data.table::last()
## x lubridate::mday()        masks data.table::mday()
## x lubridate::minute()      masks data.table::minute()
## x lubridate::month()       masks data.table::month()
## x plyr::mutate()           masks ggpubr::mutate(), dplyr::mutate()
## x tidyr::pack()            masks Matrix::pack()
## x lubridate::quarter()     masks data.table::quarter()
## x dplyr::recode()          masks car::recode()
## x plyr::rename()           masks dplyr::rename()
## x lubridate::second()      masks data.table::second()
## x lubridate::setdiff()     masks base::setdiff()
## x purrr::some()            masks car::some()
## x plyr::summarise()        masks dplyr::summarise()
## x plyr::summarize()        masks dplyr::summarize()
## x purrr::transpose()       masks data.table::transpose()
## x lubridate::union()       masks base::union()
## x tidyr::unpack()          masks Matrix::unpack()
## x lubridate::wday()        masks data.table::wday()
## x lubridate::week()        masks data.table::week()
## x lubridate::yday()        masks data.table::yday()
## x lubridate::year()        masks data.table::year()
library(vegan)
library(xlsx)
library(yaml)
###Bacteria###
unzip(zipfile = "table_filtered.16s.qza")
hdf5_biom.16s <- read_hdf5_biom("/Users/martinylab/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/1a5aa860-8bb9-46e9-8a69-2499a6e1cc21/data/feature-table.biom") 
write_biom(hdf5_biom.16s, "formatted_biom")
json_biom.16s <- read_biom("formatted_biom")
otu.table.16s <- as.data.frame(as.matrix(biom_data(json_biom.16s)))
write.table(otu.table.16s, file = "otu.table.filtered.16s.tsv", quote=FALSE, sep='\t', col.names = NA)

###Fungi###
hdf5_biom.ITS <- read_hdf5_biom("~/Google_drive/my_projects/lrgce/lrgce_seq_v2/ITS/forward_final/exported/feature-table.biom")
write_biom(hdf5_biom.ITS, "formatted_biom")
json_biom.ITS <- read_biom("formatted_biom")
otu.table.ITS <- as.data.frame(as.matrix(biom_data(json_biom.ITS)))
write.table(otu.table.ITS, file = "otu.table.filtered.ITS.tsv", quote=FALSE, sep='\t', col.names = NA)

###Plant###
plant.2015.fg <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/plant.2015.fg.xlsx")
otu.plant.2015 <- as.data.frame(plant.2015.fg, col.names    =   TRUE,   row.names   =   TRUE)
otu.plant.2015 [1:5 , 1:6]
##     sample_id      Bare    Litter AVspp      BRDI BRHO
## 1 G05LRN_2015  7.843137  4.901961     0 0.0000000    0
## 2 G05RRX_2015 11.764706 12.745098     0 0.0000000    0
## 3 G07LRX_2015  4.901961  2.941176     0 0.0000000    0
## 4 G07RRN_2015 15.686275  5.882353     0 6.8627451    0
## 5 G14LRN_2015 13.725490 14.705882     0 0.9803922    0
otu.plant.2015.1 <- otu.plant.2015[-1]
row.names(otu.plant.2015.1) <- otu.plant.2015$sample_id
otu.plant.2015.1[1:5 , 1:4]
##                  Bare    Litter AVspp      BRDI
## G05LRN_2015  7.843137  4.901961     0 0.0000000
## G05RRX_2015 11.764706 12.745098     0 0.0000000
## G07LRX_2015  4.901961  2.941176     0 0.0000000
## G07RRN_2015 15.686275  5.882353     0 6.8627451
## G14LRN_2015 13.725490 14.705882     0 0.9803922
combo.16s <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/combo.16s.xlsx")
metadata.16s <- as.data.frame(combo.16s, col.names  =   TRUE,   row.names   =   TRUE)
metadata.16s[1:5 , 1:6]
##            Sample_id vegtype_plot_treatments Ecosystem Plot Block Precipitation
## 1 14RRXAugust201616S                  G14RRX     Grass   14     6       Reduced
## 2 18LXNAugust201616S                  G18LXN     Grass   18     6       Ambient
## 3 18RXXAugust201616S                  G18RXX     Grass   18     6       Ambient
## 4 20LRXAugust201616S                  G20LRX     Grass   20     8       Reduced
## 5 22LXXAugust201616S                  G22LXX     Grass   22     8       Ambient
metadata.16s.1 <- metadata.16s[-1]
row.names(metadata.16s.1) <- metadata.16s$Sample_id
sapply(metadata.16s.1, class)
##             vegtype_plot_treatments                           Ecosystem 
##                         "character"                         "character" 
##                                Plot                               Block 
##                           "numeric"                           "numeric" 
##                       Precipitation                            Nitrogen 
##                         "character"                         "character" 
##                             eco.nit                            eco.prec 
##                         "character"                         "character" 
##                     Collection_date                              season 
##                         "character"                         "character" 
##                       eco.prec.seas                        Rainfall(in) 
##                         "character"                           "numeric" 
##                        Rainfall(cm) Rainfall_corrected_by_treatment(cm) 
##                           "numeric"                           "numeric" 
##                             shannon                            richness 
##                           "numeric"                           "numeric" 
##                            pielou_e                        average_sobs 
##                           "numeric"                           "numeric"
metadata.16s.1[1:5 , 1:6]
##                    vegtype_plot_treatments Ecosystem Plot Block Precipitation
## 14RRXAugust201616S                  G14RRX     Grass   14     6       Reduced
## 18LXNAugust201616S                  G18LXN     Grass   18     6       Ambient
## 18RXXAugust201616S                  G18RXX     Grass   18     6       Ambient
## 20LRXAugust201616S                  G20LRX     Grass   20     8       Reduced
## 22LXXAugust201616S                  G22LXX     Grass   22     8       Ambient
##                    Nitrogen
## 14RRXAugust201616S  Ambient
## 18LXNAugust201616S    Added
## 18RXXAugust201616S  Ambient
## 20LRXAugust201616S  Ambient
## 22LXXAugust201616S  Ambient
combo.ITS <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/combo.ITS.xlsx")
metadata.ITS <- as.data.frame(combo.ITS, col.names  =   TRUE,   row.names   =   TRUE)
metadata.ITS[1:5 , 1:6]
##               sample vegtype_plot_treatments Plot Block Ecosystem Precipitation
## 1 14RRXAugust2016ITS                  G14RRX   14     6     Grass       Reduced
## 2 18LXNAugust2016ITS                  G18LXN   18     6     Grass       Ambient
## 3 18RXXAugust2016ITS                  G18RXX   18     6     Grass       Ambient
## 4 20LRXAugust2016ITS                  G20LRX   20     8     Grass       Reduced
## 5 22LXXAugust2016ITS                  G22LXX   22     8     Grass       Ambient
metadata.ITS.1 <- metadata.ITS[-1]
row.names(metadata.ITS.1) <- metadata.ITS$sample
sapply(metadata.ITS.1, class)
##             vegtype_plot_treatments                                Plot 
##                         "character"                           "numeric" 
##                               Block                           Ecosystem 
##                           "numeric"                         "character" 
##                       Precipitation                            Nitrogen 
##                         "character"                         "character" 
##                             eco.nit                            eco.prec 
##                         "character"                         "character" 
##                     Collection_date                              season 
##                         "character"                         "character" 
##                        Rainfall(in)                        Rainfall(cm) 
##                           "numeric"                           "numeric" 
## Rainfall_corrected_by_treatment(cm)                             shannon 
##                           "numeric"                           "numeric" 
##                            richness                            pielou_e 
##                           "numeric"                           "numeric" 
##                        average_sobs 
##                           "numeric"
metadata.ITS.1[1:5 , 1:6]
##                    vegtype_plot_treatments Plot Block Ecosystem Precipitation
## 14RRXAugust2016ITS                  G14RRX   14     6     Grass       Reduced
## 18LXNAugust2016ITS                  G18LXN   18     6     Grass       Ambient
## 18RXXAugust2016ITS                  G18RXX   18     6     Grass       Ambient
## 20LRXAugust2016ITS                  G20LRX   20     8     Grass       Reduced
## 22LXXAugust2016ITS                  G22LXX   22     8     Grass       Ambient
##                    Nitrogen
## 14RRXAugust2016ITS  Ambient
## 18LXNAugust2016ITS    Added
## 18RXXAugust2016ITS  Ambient
## 20LRXAugust2016ITS  Ambient
## 22LXXAugust2016ITS  Ambient
plant.combo.meta.2015 <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/plant.combo.meta_2015.xlsx")
metadata.plant.2015 <- as.data.frame(plant.combo.meta.2015, col.names   =   TRUE,   row.names   =   TRUE)
metadata.plant.2015[1:5 , 1:6]
##     sample_id   plot block year ecosystem nitrogen
## 1 G05LRN_2015 G05LRN     2 2015     grass    added
## 2 G05RRX_2015 G05RRX     2 2015     grass  ambient
## 3 G07LRX_2015 G07LRX     3 2015     grass  ambient
## 4 G07RRN_2015 G07RRN     3 2015     grass    added
## 5 G14LRN_2015 G14LRN     6 2015     grass    added
metadata.plant.2015.1 <- metadata.plant.2015[-1]
row.names(metadata.plant.2015.1) <- metadata.plant.2015$sample_id
sapply(metadata.plant.2015.1, class)
##                plot               block                year           ecosystem 
##         "character"           "numeric"           "numeric"         "character" 
##            nitrogen       precipitation        rainfall(cm)            eco.prec 
##         "character"         "character"           "numeric"         "character" 
##             eco.nit  NonNativeForbCover NonNativeGrassCover     NativeForbCover 
##         "character"           "numeric"           "numeric"           "numeric" 
##    NativeGrassCover    NativeShrubCover                Bare              Litter 
##           "numeric"           "numeric"           "numeric"           "numeric" 
##   total_grass_cover TotalNonNativeCover    TotalNativeCover 
##           "numeric"           "numeric"           "numeric"
metadata.plant.2015.1[1:5 , 1:6]
##               plot block year ecosystem nitrogen precipitation
## G05LRN_2015 G05LRN     2 2015     grass    added       reduced
## G05RRX_2015 G05RRX     2 2015     grass  ambient       reduced
## G07LRX_2015 G07LRX     3 2015     grass  ambient       reduced
## G07RRN_2015 G07RRN     3 2015     grass    added       reduced
## G14LRN_2015 G14LRN     6 2015     grass    added       reduced
taxonomy.16s <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/taxonomy_16s.xlsx")
taxaid.16s <- as.data.frame(taxonomy.16s, col.names =   TRUE,   row.names   =   TRUE)
taxaid.16s[1:2 , 1:2]
##                              OTUID
## 1 9ccef5c0259d44ee21a37f1838ef585b
## 2 6f301efd81a9f8bfb9c3d39ff3700c8a
##                                                                                                                           taxonomy
## 1 k__Bacteria; p__Actinobacteria; c__Actinobacteria; o__Actinomycetales; f__Streptomycetaceae; g__Streptomyces; s__reticuliscabiei
## 2                                         k__Bacteria; p__Bacteroidetes; c__Sphingobacteriia; o__Sphingobacteriales; f__; g__; s__
taxaid.16s.1 <-taxaid.16s[-1]
row.names(taxaid.16s.1) <- taxaid.16s$OTUID
sapply(taxaid.16s.1, class)
##    taxonomy  confidence 
## "character"   "numeric"
taxaid.16s.1[1:2 , 1:2]
##                                                                                                                                                          taxonomy
## 9ccef5c0259d44ee21a37f1838ef585b k__Bacteria; p__Actinobacteria; c__Actinobacteria; o__Actinomycetales; f__Streptomycetaceae; g__Streptomyces; s__reticuliscabiei
## 6f301efd81a9f8bfb9c3d39ff3700c8a                                         k__Bacteria; p__Bacteroidetes; c__Sphingobacteriia; o__Sphingobacteriales; f__; g__; s__
##                                  confidence
## 9ccef5c0259d44ee21a37f1838ef585b  0.9828906
## 6f301efd81a9f8bfb9c3d39ff3700c8a  0.9999996
taxonomy.ITS <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/taxonomy_ITS.xlsx")
taxaid.ITS <- as.data.frame(taxonomy.ITS, col.names =   TRUE,   row.names   =   TRUE)
taxaid.ITS[1:2 , 1:2]
##                             otu_id
## 1 7e80c7e5f5813d51bcb914eae7555273
## 2 1094975b0bdad9f93c92f4b3b2e6df15
##                                                                                                                taxon
## 1                                                                                             k__Fungi;p__Ascomycota
## 2 k__Fungi;p__Ascomycota;c__Leotiomycetes;o__Helotiales;f__Helotiaceae;g__Articulospora;s__Articulospora_proliferata
taxaid.ITS.1 <-taxaid.ITS[-1]
row.names(taxaid.ITS.1) <- taxaid.ITS$otu_id
sapply(taxaid.ITS.1, class)
##       taxon  confidence 
## "character"   "numeric"
taxaid.ITS.1[1:2 , 1:2]
##                                                                                                                                               taxon
## 7e80c7e5f5813d51bcb914eae7555273                                                                                             k__Fungi;p__Ascomycota
## 1094975b0bdad9f93c92f4b3b2e6df15 k__Fungi;p__Ascomycota;c__Leotiomycetes;o__Helotiales;f__Helotiaceae;g__Articulospora;s__Articulospora_proliferata
##                                  confidence
## 7e80c7e5f5813d51bcb914eae7555273  0.8120831
## 1094975b0bdad9f93c92f4b3b2e6df15  0.9517936
###Bacteria###
otu.plus.taxaid.16s <- as.data.frame(merge(taxaid.16s.1, otu.table.16s, by.x = "row.names", by.y = "row.names"))
otu.plus.taxaid.filtered.16s <- otu.plus.taxaid.16s[!grepl("chloroplast", otu.plus.taxaid.16s$`taxon`),]
otu.plus.taxaid.filtered.16s <- otu.plus.taxaid.filtered.16s[!grepl("mitochondria",otu.plus.taxaid.filtered.16s$`taxon`),]
write.table(otu.plus.taxaid.filtered.16s, file = "otu.taxaid.filtered.16s.tsv", quote=FALSE, sep='\t', col.names = NA)
row.names(otu.plus.taxaid.filtered.16s) <- otu.plus.taxaid.filtered.16s$Row.names
otu.clean.16s <- as.data.frame(t(otu.plus.taxaid.filtered.16s[,4:ncol(otu.plus.taxaid.filtered.16s)]))
write.table(otu.clean.16s, file = "otu.clean.16s.tsv", quote=FALSE, sep='\t', col.names = NA)

read.depth.barplot.16s <- barplot(sort(rowSums(otu.clean.16s)), ylim = c(0, max(rowSums(otu.clean.16s))), 
                          xlim = c(0, NROW(otu.clean.16s)), col = "Blue", ylab = "Read Depth", xlab = "Sample") 

unrarfied.curves.16s <-rarecurve(otu.clean.16s, step = 1000, label = FALSE) 

set.seed(seed = 999)                  
rared.otu.16s <- as.data.frame((rrarefy.perm(otu.clean.16s, sample = 1090, n = 300, round.out = T)))#because rrarefy.perm randomly picks OTUs at a specified sequencing depth, the total number of OTUs may change with each run of the function.
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  2  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  3  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  4  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  5  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  6  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  7  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  8  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  9  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  10  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  11  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  12  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  13  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  14  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  15  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  16  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  17  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  18  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  19  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  20  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  21  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  22  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  23  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  24  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  25  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  26  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  27  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  28  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  29  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  30  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  31  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  32  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  33  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  34  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  35  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  36  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  37  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  38  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  39  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  40  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  41  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  42  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  43  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  44  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  45  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  46  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  47  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  48  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  49  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  50  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  51  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  52  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  53  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  54  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  55  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  56  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  57  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  58  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  59  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  60  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  61  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  62  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  63  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  64  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  65  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  66  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  67  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  68  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  69  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  70  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  71  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  72  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  73  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  74  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  75  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  76  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  77  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  78  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  79  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  80  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  81  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  82  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  83  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  84  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  85  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  86  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  87  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  88  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  89  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  90  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  91  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  92  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  93  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  94  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  95  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  96  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  97  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  98  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  99  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  100  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  101  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  102  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  103  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  104  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  105  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  106  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  107  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  108  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  109  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  110  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  111  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  112  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  113  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  114  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  115  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  116  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  117  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  118  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  119  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  120  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  121  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  122  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  123  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  124  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  125  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  126  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  127  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  128  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  129  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  130  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  131  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  132  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  133  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  134  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  135  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  136  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  137  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  138  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  139  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  140  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  141  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  142  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  143  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  144  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  145  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  146  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  147  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  148  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  149  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  150  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  151  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  152  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  153  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  154  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  155  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  156  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  157  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  158  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  159  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  160  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  161  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  162  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  163  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  164  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  165  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  166  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  167  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  168  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  169  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  170  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  171  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  172  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  173  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  174  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  175  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  176  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  177  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  178  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  179  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  180  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  181  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  182  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  183  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  184  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  185  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  186  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  187  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  188  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  189  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  190  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  191  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  192  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  193  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  194  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  195  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  196  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  197  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  198  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  199  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  200  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  201  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  202  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  203  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  204  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  205  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  206  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  207  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  208  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  209  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  210  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  211  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  212  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  213  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  214  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  215  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  216  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  217  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  218  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  219  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  220  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  221  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  222  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  223  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  224  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  225  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  226  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  227  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  228  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  229  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  230  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  231  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  232  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  233  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  234  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  235  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  236  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  237  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  238  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  239  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  240  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  241  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  242  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  243  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  244  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  245  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  246  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  247  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  248  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  249  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  250  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  251  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  252  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  253  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  254  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  255  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  256  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  257  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  258  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  259  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  260  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  261  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  262  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  263  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  264  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  265  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  266  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  267  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  268  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  269  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  270  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  271  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  272  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  273  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  274  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  275  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  276  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  277  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  278  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  279  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  280  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  281  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  282  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  283  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  284  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  285  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  286  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  287  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  288  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  289  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  290  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  291  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  292  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  293  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  294  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  295  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  296  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  297  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  298  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  299  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  300  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
warnings(rared.otu.16s)
rared.otu.16s <- as.data.frame(rared.otu.16s[rowSums(rared.otu.16s) >= 1090-(1090*.1), colSums(rared.otu.16s) >= 1])

rarefied.curves.16s <- rarecurve(rared.otu.16s, step = 300)

rared.otu.16s.t <- t(rared.otu.16s)
rared.otu.taxaid.16s <- merge(rared.otu.16s.t, taxaid.16s.1, by.x = "row.names", by.y = "row.names")
write.table(rared.otu.taxaid.16s, file = "rared.otu.taxaid.16s.tsv", quote=FALSE, sep='\t', col.names = NA)

###Fungi###
otu.plus.taxaid.ITS <- as.data.frame(merge(taxaid.ITS.1, otu.table.ITS, by.x = "row.names", by.y = "row.names"))
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.ITS[!grepl("k__Chromista", otu.plus.taxaid.ITS$`taxon`),]
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.filtered.ITS[!grepl("k__Rhizaria", otu.plus.taxaid.filtered.ITS$`taxon`),]
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.filtered.ITS[!grepl("k__Plantae", otu.plus.taxaid.filtered.ITS$`taxon`),]
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.filtered.ITS[!grepl("k__Alveolata", otu.plus.taxaid.filtered.ITS$`taxon`),]
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.filtered.ITS[!grepl("k__Metazoa", otu.plus.taxaid.filtered.ITS$`taxon`),]
otu.plus.taxaid.filtered.ITS <- otu.plus.taxaid.filtered.ITS[!grepl("k__Fungi;p__unidentified", otu.plus.taxaid.filtered.ITS$`taxon`),]
write.table(otu.plus.taxaid.filtered.ITS, file = "otu.taxaid.filtered.ITS.tsv", quote=FALSE, sep='\t', col.names = NA)
row.names(otu.plus.taxaid.filtered.ITS) <- otu.plus.taxaid.filtered.ITS$Row.names
otu.clean.ITS <- as.data.frame(t(otu.plus.taxaid.filtered.ITS[,4:ncol(otu.plus.taxaid.filtered.ITS)]))
write.table(otu.clean.ITS, file = "otu.clean.ITS.tsv", quote=FALSE, sep='\t', col.names = NA)

read.depth.barplot.ITS <- barplot(sort(rowSums(otu.clean.ITS)), ylim = c(0, max(rowSums(otu.clean.ITS))), 
                          xlim = c(0, NROW(otu.clean.ITS)), col = "Blue", ylab = "Read Depth", xlab = "Sample") 

unrarfied.curves.ITS <-rarecurve(otu.clean.ITS, step = 1000, label = FALSE) 

set.seed(seed = 999)                  
rared.otu.ITS <- as.data.frame((rrarefy.perm(otu.clean.ITS, sample = 1064, n = 300, round.out = T)))
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  2  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  3  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  4  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  5  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  6  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  7  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  8  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  9  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  10  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  11  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  12  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  13  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  14  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  15  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  16  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  17  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  18  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  19  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  20  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  21  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  22  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  23  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  24  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  25  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  26  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  27  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  28  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  29  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  30  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  31  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  32  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  33  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  34  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  35  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  36  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  37  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  38  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  39  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  40  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  41  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  42  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  43  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  44  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  45  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  46  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  47  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  48  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  49  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  50  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  51  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  52  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  53  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  54  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  55  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  56  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  57  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  58  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  59  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  60  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  61  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  62  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  63  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  64  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  65  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  66  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  67  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  68  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  69  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  70  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  71  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  72  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  73  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  74  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  75  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  76  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  77  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  78  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  79  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  80  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  81  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  82  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  83  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  84  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  85  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  86  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  87  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  88  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  89  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  90  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  91  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  92  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  93  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  94  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  95  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  96  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  97  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  98  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  99  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  100  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  101  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  102  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  103  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  104  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  105  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  106  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  107  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  108  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  109  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  110  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  111  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  112  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  113  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  114  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  115  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  116  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  117  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  118  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  119  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  120  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  121  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  122  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  123  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  124  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  125  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  126  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  127  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  128  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  129  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  130  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  131  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  132  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  133  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  134  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  135  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  136  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  137  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  138  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  139  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  140  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  141  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  142  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  143  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  144  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  145  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  146  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  147  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  148  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  149  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  150  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  151  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  152  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  153  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  154  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  155  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  156  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  157  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  158  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  159  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  160  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  161  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  162  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  163  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  164  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  165  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  166  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  167  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  168  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  169  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  170  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  171  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  172  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  173  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  174  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  175  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  176  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  177  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  178  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  179  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  180  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  181  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  182  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  183  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  184  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  185  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  186  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  187  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  188  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  189  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  190  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  191  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  192  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  193  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  194  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  195  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  196  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  197  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  198  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  199  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  200  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  201  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  202  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  203  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  204  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  205  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  206  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  207  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  208  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  209  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  210  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  211  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  212  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  213  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  214  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  215  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  216  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  217  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  218  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  219  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  220  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  221  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  222  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  223  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  224  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  225  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  226  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  227  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  228  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  229  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  230  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  231  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  232  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  233  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  234  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  235  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  236  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  237  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  238  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  239  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  240  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  241  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  242  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  243  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  244  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  245  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  246  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  247  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  248  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  249  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  250  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  251  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  252  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  253  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  254  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  255  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  256  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  257  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  258  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  259  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  260  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  261  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  262  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  263  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  264  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  265  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  266  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  267  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  268  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  269  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  270  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  271  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  272  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  273  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  274  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  275  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  276  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  277  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  278  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  279  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  280  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  281  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  282  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  283  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  284  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  285  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  286  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  287  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  288  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  289  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  290  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  291  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  292  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  293  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  294  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  295  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  296  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  297  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  298  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  299  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
## Permutation  300  out of  300
## Warning in rrarefy(x, sample): some row sums < 'sample' and are not rarefied
warnings(rared.otu.ITS)
rared.otu.ITS <- as.data.frame(rared.otu.ITS[rowSums(rared.otu.ITS) >= 1064-(1064*.1), colSums(rared.otu.ITS) >= 1])

rarefied.curves.ITS <- rarecurve(rared.otu.ITS, step = 300)

rared.otu.ITS.t <- t(rared.otu.ITS)
rared.otu.taxaid.ITS <- merge(rared.otu.ITS.t, taxaid.ITS.1, by.x = "row.names", by.y = "row.names")
write.table(rared.otu.taxaid.ITS, file = "rared.otu.taxaid.ITS.tsv", quote=FALSE, sep='\t', col.names = NA)
#Create function to read qza file into R
read_qza <- function(file, tmp, rm) {
  
  if(missing(tmp)){tmp <- tempdir()}
  if(missing(file)){stop("Path to artifact (.qza) not provided")}
  if(missing(rm)){rm=TRUE} #remove the decompressed object from tmp
  
  unzip(file, exdir=tmp)
  unpacked<-unzip(file, exdir=tmp, list=TRUE)
  
  artifact<-read_yaml(paste0(tmp,"/", paste0(gsub("/..+","", unpacked$Name[1]),"/metadata.yaml"))) #start by loading in the metadata not assuming it will be first file listed
  artifact$contents<-data.frame(files=unpacked)
  artifact$contents$size=sapply(paste0(tmp, "/", artifact$contents$files), file.size)
  artifact$version=read.table(paste0(tmp,"/",artifact$uuid, "/VERSION"))
  
#if(sum(artifact$version$V2==c("2","4","2018.4.0"))!=3){warning("Artifact was not generated with Qiime2 2018.4, if data is not successfully imported, please report here github.com/jbisanz/qiime2R/issues")}#check version and throw warning if new format
  
#get data dependent on format
  if(grepl("BIOMV", artifact$format)){
    suppressWarnings(artifact$data<-as(biom_data(read_biom(paste0(tmp, "/", artifact$uui,"/data/feature-table.biom"))),"matrix")) #suppressing warning about \n characters
  } else if (artifact$format=="NewickDirectoryFormat"){
    artifact$data<-read.tree(paste0(tmp,"/",artifact$uuid,"/data/tree.nwk"))
  } else if (artifact$format=="DistanceMatrixDirectoryFormat") {
    artifact$data<-as.dist(read.table(paste0(tmp,"/", artifact$uuid, "/data/distance-matrix.tsv"), header=TRUE, row.names=1))
  } else if (grepl("StatsDirFmt", artifact$format)) {
    if(paste0(artifact$uuid, "/data/stats.csv") %in% artifact$contents$files.Name){artifact$data<-read.csv(paste0(tmp,"/", artifact$uuid, "/data/stats.csv"), header=TRUE, row.names=1)}
    if(paste0(artifact$uuid, "/data/stats.tsv") %in% artifact$contents$files.Name){artifact$data<-read.table(paste0(tmp,"/", artifact$uuid, "/data/stats.tsv"), header=TRUE, row.names=1, sep='\t')} #can be tsv or csv
  } else if (artifact$format=="TSVTaxonomyDirectoryFormat"){
    artifact$data<-read.table(paste0(tmp,"/", artifact$uuid, "/data/taxonomy.tsv"), sep='\t', header=TRUE, quote="")
  } else if (artifact$format=="OrdinationDirectoryFormat"){
    
    linesplit<-suppressWarnings(readLines(paste0(tmp,"/", artifact$uuid, "/data/ordination.txt")))
    linesplit<-linesplit[sapply(linesplit, function(x) x!="")]
    
    for (i in 1:length(linesplit)){
      if(grepl("^Eigvals\\t|^Proportion explained\\t|^Species\\t|^Site\\t|^Biplot\\t|^Site constraints\\t", linesplit[i])){
        curfile=strsplit(linesplit[i],"\t")[[1]][1]
      } else {
        write(linesplit[i], paste0(tmp,"/", artifact$uuid, "/data/",curfile,".tmp"), append=TRUE)
      }
    }
    
    for (outs in list.files(paste0(tmp,"/", artifact$uuid,"/data"), full.names = TRUE, pattern = "\\.tmp")){
      NewLab<-gsub(" ", "", toTitleCase(gsub("\\.tmp", "", basename(outs))))
      artifact$data[[NewLab]]<-read.table(outs,sep='\t', header=FALSE)
      if(NewLab %in% c("Eigvals","ProportionExplained")){colnames(artifact$data[[NewLab]])<-paste0("PC",1:ncol(artifact$data[[NewLab]]))}
      if(NewLab %in% c("Site","SiteConstraints")){colnames(artifact$data[[NewLab]])<-c("SampleID", paste0("PC",1:(ncol(artifact$data[[NewLab]])-1)))}
      if(NewLab %in% c("Species")){colnames(artifact$data[[NewLab]])<-c("FeatureID", paste0("PC",1:(ncol(artifact$data[[NewLab]])-1)))}
    }
    
    artifact$data$Vectors<-artifact$data$Site #Rename Site to Vectors so this matches up with the syntax used in the tutorials
    artifact$data$Site<-NULL
    
  } else if (artifact$format=="DNASequencesDirectoryFormat") {
    artifact$data<-readDNAStringSet(paste0(tmp,"/",artifact$uuid,"/data/dna-sequences.fasta"))
  } else if (artifact$format=="AlignedDNASequencesDirectoryFormat") {
    artifact$data<-readDNAMultipleAlignment(paste0(tmp,"/",artifact$uuid,"/data/aligned-dna-sequences.fasta"))
  } else if (grepl("EMPPairedEndDirFmt|EMPSingleEndDirFmt|FastqGzFormat|MultiplexedPairedEndBarcodeInSequenceDirFmt|MultiplexedSingleEndBarcodeInSequenceDirFmt|PairedDNASequencesDirectoryFormat|SingleLanePerSamplePairedEndFastqDirFmt|SingleLanePerSampleSingleEndFastqDirFmt", artifact$format)) {
    artifact$data<-data.frame(files=list.files(paste0(tmp,"/", artifact$uuid,"/data")))
    artifact$data$size<-format(sapply(artifact$data$files, function(x){file.size(paste0(tmp,"/",artifact$uuid,"/data/",x))}, simplify = TRUE))
  } else if (artifact$format=="AlphaDiversityDirectoryFormat") {
    artifact$data<-read.table(paste0(tmp, "/", artifact$uuid, "/data/alpha-diversity.tsv"))
  } else {
    message("Format not supported, only a list of internal files and provenance is being imported.")
    artifact$data<-list.files(paste0(tmp,"/",artifact$uuid, "/data"))
  }
  
  pfiles<-paste0(tmp,"/", grep("..+provenance/..+action.yaml", unpacked$Name, value=TRUE))
  artifact$provenance<-lapply(pfiles, read_yaml)
  names(artifact$provenance)<-grep("..+provenance/..+action.yaml", unpacked$Name, value=TRUE)
  if(rm==TRUE){unlink(paste0(tmp,"/", artifact$uuid), recursive=TRUE)}
  return(artifact)
}

###Bacteria###
taxonomy16s<-read_qza("taxonomy16s.qza", "./tmp", rm=T)
taxonomy16s$uuid
## [1] "67172321-d90d-436c-8328-b24d8743c424"
taxtable16s<-taxonomy16s$data %>% as.tibble() %>% separate(Taxon, sep="; ", c("Kingdom","Phylum","Class","Order","Family","Genus","Species")) #convert the table into a tabular split version
## Warning: `as.tibble()` is deprecated as of tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: Expected 7 pieces. Missing pieces filled with `NA` in 749 rows [3, 10,
## 12, 15, 16, 19, 24, 29, 34, 37, 40, 42, 45, 53, 55, 57, 60, 61, 62, 63, ...].
taxtable16s
## # A tibble: 1,866 x 9
##    Feature.ID   Kingdom  Phylum  Class  Order  Family  Genus  Species Confidence
##    <fct>        <chr>    <chr>   <chr>  <chr>  <chr>   <chr>  <chr>        <dbl>
##  1 9ccef5c0259… k__Bact… p__Act… c__Ac… o__Ac… f__Str… g__St… s__ret…      0.983
##  2 6f301efd81a… k__Bact… p__Bac… c__Sp… o__Sp… f__     g__    s__          1.00 
##  3 9b1857d7d74… k__Bact… <NA>    <NA>   <NA>   <NA>    <NA>   <NA>         0.960
##  4 5888bd387cc… k__Bact… p__Pro… c__De… o__My… f__     g__    s__          0.798
##  5 f26c82353ae… k__Bact… p__Aci… c__[C… o__RB… f__Ell… g__    s__          1.00 
##  6 8e9ec770809… k__Bact… p__Bac… c__Fl… o__Fl… f__[We… g__Ch… s__          1.00 
##  7 6419cb3e5b3… k__Bact… p__Act… c__Ac… o__Ac… f__San… g__Sa… s__          0.993
##  8 fb33e463293… k__Bact… p__Pro… c__Al… o__Sp… f__Sph… g__Sp… s__          0.907
##  9 cadb3361e0a… k__Bact… p__Act… c__Ac… o__Ac… f__Kin… g__Ki… s__          0.998
## 10 aa97c16a0a1… k__Bact… p__Pro… c__Be… o__Bu… f__Com… g__Li… <NA>         0.703
## # … with 1,856 more rows
taxtable16s <- as.data.frame(taxtable16s)

#remove p__ as a precursor to phylum name 
taxtable16s$Phylum <- gsub("p__", "", taxtable16s$Phylum)
taxtable16s$Class <- gsub("c__", "", taxtable16s$Class)
taxtable16s$Order <- gsub("o__", "", taxtable16s$Order)
taxtable16s$Family <- gsub("f__", "", taxtable16s$Family)
taxtable16s$Genus <- gsub("g__", "", taxtable16s$Genus)
taxtable16s$Species <- gsub("s__", "", taxtable16s$Species)

#this will turn all empty cell into NA
taxtable16s <- taxtable16s %>%
  mutate_at(vars(colnames(.)),
            .funs = funs(ifelse(.=="", NA, as.character(.))))
## Warning: funs() is soft deprecated as of dplyr 0.8.0
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once per session.
#taxtable was imported using qiime2R, with already separated levels
row.names(taxtable16s) <- taxtable16s$Feature.ID #turn feature IDs into row.names
taxtable16s$Feature.ID <- NULL #remove column "Feature ID"

#merge OTU table with taxonomy table
otu.raref.taxonomy.16s <- as.data.frame(merge(taxtable16s, rared.otu.16s.t, by.x = "row.names", by.y = "row.names"))
row.names(otu.raref.taxonomy.16s) <- otu.raref.taxonomy.16s$Row.names #move column "Row.names" (first column) to row.names this was created during the merge
otu.raref.taxonomy.16s$Row.names <- NULL #remove column "Row.names"
otu.raref.taxonomy.16s$Confidence <- NULL # remove column "Confidence"
otu.raref.taxonomy.16s.ss <- subset(otu.raref.taxonomy.16s, select=c(8:195,1:7)) # move Taxa information to the end of all columns
write.table(otu.raref.taxonomy.16s.ss, "otu_raref_taxonomy.16s.csv", quote=FALSE, sep=',') #use this table to calculate top number of ESVs  and their sequence contributions.
colSums(otu.raref.taxonomy.16s.ss[, c(1:188)]) # check if samples are RAREFIED!!
##  14LRNDecember201616S      14LRNJune201716S     14LRNMarch201716S 
##                  1090                  1090                  1097 
##     14LRNMarch201816S    14RRXAugust201616S  14RRXDecember201616S 
##                  1090                  1088                  1088 
##  14RRXDecember201716S      14RRXJune201716S     14RRXMarch201716S 
##                  1091                  1090                  1092 
##     14RRXMarch201816S 14RRXSeptember201716S    18LXNAugust201616S 
##                  1093                  1089                  1087 
##  18LXNDecember201616S  18LXNDecember201716S      18LXNJune201716S 
##                  1093                  1089                  1090 
##     18LXNMarch201716S     18LXNMarch201816S 18LXNSeptember201716S 
##                  1091                  1089                  1091 
##    18RXXAugust201616S  18RXXDecember201716S      18RXXJune201716S 
##                  1090                  1091                  1089 
##     18RXXMarch201716S     18RXXMarch201816S 18RXXSeptember201716S 
##                  1089                  1087                  1090 
##    20LRXAugust201616S  20LRXDecember201616S  20LRXDecember201716S 
##                  1090                  1087                  1090 
##      20LRXJune201716S     20LRXMarch201716S     20LRXMarch201816S 
##                  1088                  1091                  1090 
## 20LRXSeptember201716S  20RRNDecember201616S  20RRNDecember201716S 
##                  1089                  1087                  1091 
##     20RRNMarch201716S     20RRNMarch201816S    22LXXAugust201616S 
##                  1092                  1088                  1091 
##  22LXXDecember201616S  22LXXDecember201716S     22LXXMarch201716S 
##                  1092                  1090                  1092 
##     22LXXMarch201816S    22RXNAugust201616S  22RXNDecember201616S 
##                  1089                  1088                  1091 
##  22RXNDecember201716S      22RXNJune201716S     22RXNMarch201716S 
##                  1090                  1088                  1089 
##     22RXNMarch201816S 22RXNSeptember201716S    25LRXAugust201616S 
##                  1088                  1088                  1089 
##  25LRXDecember201616S  25LRXDecember201716S     25LRXMarch201716S 
##                  1086                  1089                  1090 
##     25LRXMarch201816S 25LRXSeptember201716S    25RRNAugust201616S 
##                  1090                  1092                  1082 
##  25RRNDecember201616S  25RRNDecember201716S      25RRNJune201716S 
##                  1089                  1088                  1091 
##     25RRNMarch201716S     25RRNMarch201816S 25RRNSeptember201716S 
##                  1092                  1092                  1092 
##    27LXNAugust201616S  27LXNDecember201616S  27LXNDecember201716S 
##                  1097                  1096                  1088 
##      27LXNJune201716S     27LXNMarch201716S     27LXNMarch201816S 
##                  1091                  1089                  1091 
## 27LXNSeptember201716S    27RXXAugust201616S  27RXXDecember201616S 
##                  1090                  1086                  1090 
##  27RXXDecember201716S      27RXXJune201716S     27RXXMarch201716S 
##                  1087                  1091                  1087 
##     27RXXMarch201816S 27RXXSeptember201716S    32LXNAugust201616S 
##                  1087                  1090                  1090 
##  32LXNDecember201616S     32LXNMarch201716S     32LXNMarch201816S 
##                  1089                  1088                  1090 
## 32LXNSeptember201716S    32RXXAugust201616S  32RXXDecember201616S 
##                  1097                  1091                  1090 
##  32RXXDecember201716S      32RXXJune201716S     32RXXMarch201716S 
##                  1090                  1092                  1092 
##     32RXXMarch201816S 32RXXSeptember201716S    35LRXAugust201616S 
##                  1089                  1090                  1089 
##  35LRXDecember201616S  35LRXDecember201716S      35LRXJune201716S 
##                  1093                  1091                  1090 
##     35LRXMarch201716S     35LRXMarch201816S 35LRXSeptember201716S 
##                  1090                  1088                  1090 
##    35RRNAugust201616S  35RRNDecember201616S     35RRNMarch201716S 
##                  1090                  1091                  1090 
##     35RRNMarch201816S 35RRNSeptember201716S    45LRXAugust201616S 
##                  1093                  1091                  1089 
##  45LRXDecember201616S      45LRXJune201716S     45LRXMarch201716S 
##                  1089                  1091                  1091 
##     45LRXMarch201816S 45LRXSeptember201716S    45RRNAugust201616S 
##                  1090                  1089                  1092 
##  45RRNDecember201616S  45RRNDecember201716S      45RRNJune201716S 
##                  1093                  1089                  1088 
##     45RRNMarch201716S     45RRNMarch201816S 45RRNSeptember201716S 
##                  1090                  1092                  1091 
##    46LXNAugust201616S  46LXNDecember201616S  46LXNDecember201716S 
##                  1089                  1090                  1092 
##      46LXNJune201716S     46LXNMarch201716S     46LXNMarch201816S 
##                  1087                  1093                  1090 
##    46RXXAugust201616S  46RXXDecember201616S      46RXXJune201716S 
##                  1092                  1091                  1091 
##     46RXXMarch201716S     46RXXMarch201816S 46RXXSeptember201716S 
##                  1089                  1092                  1090 
##    47LRNAugust201616S  47LRNDecember201616S      47LRNJune201716S 
##                  1090                  1088                  1090 
##     47LRNMarch201716S     47LRNMarch201816S 47LRNSeptember201716S 
##                  1091                  1090                  1090 
##    47RRXAugust201616S  47RRXDecember201616S      47RRXJune201716S 
##                  1090                  1091                  1093 
##     47RRXMarch201716S     47RRXMarch201816S 47RRXSeptember201716S 
##                  1090                  1086                  1090 
##    48LXXAugust201616S  48LXXDecember201616S  48LXXDecember201716S 
##                  1089                  1088                  1087 
##      48LXXJune201716S     48LXXMarch201716S     48LXXMarch201816S 
##                  1088                  1089                  1088 
##    48RXNAugust201616S  48RXNDecember201616S      48RXNJune201716S 
##                  1087                  1086                  1090 
##     48RXNMarch201716S     48RXNMarch201816S 48RXNSeptember201716S 
##                  1092                  1091                  1092 
##   4LXXDecember201616S   4LXXDecember201716S      4LXXMarch201716S 
##                  1090                  1090                  1089 
##      4LXXMarch201816S     4RXNAugust201616S   4RXNDecember201616S 
##                  1089                  1088                  1089 
##   4RXNDecember201716S       4RXNJune201716S      4RXNMarch201716S 
##                  1089                  1090                  1087 
##      4RXNMarch201816S  4RXNSeptember201716S   5LRNDecember201616S 
##                  1090                  1089                  1091 
##   5LRNDecember201716S      5LRNMarch201716S      5LRNMarch201816S 
##                  1090                  1089                  1091 
##  5LRNSeptember201716S     5RRXAugust201616S   5RRXDecember201616S 
##                  1090                  1089                  1091 
##       5RRXJune201716S      5RRXMarch201716S      5RRXMarch201816S 
##                  1089                  1089                  1091 
##  5RRXSeptember201716S   7LRXDecember201616S   7LRXDecember201716S 
##                  1091                  1093                  1090 
##      7LRXMarch201716S      7LRXMarch201816S  7LRXSeptember201716S 
##                  1092                  1090                  1090 
##   7RRNDecember201616S   7RRNDecember201716S       7RRNJune201716S 
##                  1091                  1089                  1091 
##      7RRNMarch201716S      7RRNMarch201816S  7RRNSeptember201716S 
##                  1092                  1091                  1090 
##     8LXXAugust201616S   8LXXDecember201616S   8LXXDecember201716S 
##                  1091                  1089                  1091 
##      8LXXMarch201716S      8LXXMarch201816S  8LXXSeptember201716S 
##                  1090                  1088                  1090 
##     8RXNAugust201616S   8RXNDecember201616S 
##                  1088                  1090
# Check the number of unique OTUs at different levels
phylum.16s <- as.vector(unique(otu.raref.taxonomy.16s.ss$Phylum[!is.na(otu.raref.taxonomy.16s.ss$Phylum)]))
length(phylum.16s) #15 phyla
## [1] 15
class.16s <- as.vector(unique(otu.raref.taxonomy.16s.ss$Class[!is.na(otu.raref.taxonomy.16s.ss$Class)]))
length(class.16s) #36 unique bacterial classes (NA not counted) 
## [1] 34
order.16s <- as.vector(unique(otu.raref.taxonomy.16s.ss$Order[!is.na(otu.raref.taxonomy.16s.ss$Order)]))
length(order.16s) #55 unique bacterial orders (na not counted) 
## [1] 53
family.16s <- as.vector(unique(otu.raref.taxonomy.16s.ss$Family[!is.na(otu.raref.taxonomy.16s.ss$Family)]))
length(family.16s) #97 unique families (na not counted)
## [1] 96
genus.16s <- as.vector(unique(otu.raref.taxonomy.16s.ss$Genus[!is.na(otu.raref.taxonomy.16s.ss$Genus)]))
length(genus.16s) #142 unique bacterial genus' (na not counted) 
## [1] 141
#aggregate taxonomical IDs on a certain level by sample (here at Order level)
genus.16s.1 <- aggregate(. ~Genus, data = otu.raref.taxonomy.16s.ss[c(1:(length(names(otu.raref.taxonomy.16s.ss)) - 7), match("Genus", names(otu.raref.taxonomy.16s.ss)))], FUN = sum, na.omit=TRUE)

#Order column is turned into rownames and the empty cell is called "unidentified" (o__ gsub into "")
rownames(genus.16s.1) <- genus.16s.1$Genus
genus.16s.1$Genus <- NULL

#turn table into dataframe and transpose
genus.16s.1 <- as.data.frame(t(genus.16s.1))
write.xlsx(genus.16s.1, file = "genera.16s.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names=TRUE)

# Merge a metadata column here (I just made my own metatdata column)
matching.metadata.gen.16s <- metadata.16s.1[c(match(sort(row.names(genus.16s.1)), sort(row.names(metadata.16s.1)))),]
merged.meta.gen.16s <- merge(genus.16s.1, matching.metadata.gen.16s, by.x = "row.names", by.y = "row.names")
colnames(merged.meta.gen.16s)[colnames(merged.meta.gen.16s)=="Row.names"] <- "SampleID"
merged.meta.gen.16s <- as.data.frame(merged.meta.gen.16s)
rownames(merged.meta.gen.16s) <- merged.meta.gen.16s$SampleID
merged.meta.gen.16s$SampleID <- NULL

#aggregate only the X number of the level of taxonomy chosen (from the sequence summaries above) with the metadata selected in this case plot Treatment
genera.agg.16s <- aggregate(merged.meta.gen.16s[, 1:142], list(merged.meta.gen.16s$eco.prec), mean) #genus
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
genera.agg.16s <- as.data.frame(t(genera.agg.16s)) # transpose table to have variables as columns
names(genera.agg.16s) <- as.character(unlist(genera.agg.16s["Group.1", ])) # transposing causes the treatments named as "Group.1", remove this
genera.agg.16s <- genera.agg.16s[-1,] #final removal of "Group.1" 
for(i in c(1:ncol(genera.agg.16s))) {
  genera.agg.16s[,i] <- as.numeric(as.character(genera.agg.16s[,i]))
}
write.xlsx(genera.agg.16s, file = "plot.genera.16s.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names=TRUE)

#calculate the number of sequences per level (change the divisor according to above calculation for unique OTUs per level)
totalseqs.16s <- colSums(genera.agg.16s) 
avgnumseqs.16s <- sum(totalseqs.16s)/142 #142 is the total number of unique genera with NA removed

# CHECK FOR THE MOST ABUNDANT of classifaction level selected HERE!!!!
genera.agg.16s$Genus <- row.names(genera.agg.16s)
genera.agg.16s$Genus <- factor(genera.agg.16s$Genus)
sumsum.16s <- colSums(genera.agg.16s[1:(length(names(genera.agg.16s))-1)])
samp.16s <- names(genera.agg.16s)[1:(length(names(genera.agg.16s))-1)]
genera.agg.long.16s <- melt(genera.agg.16s, id.vars = "Genus", measure.vars = c("CSS_Ambient","CSS_Reduced", "Grass_Ambient","Grass_Reduced"))
## Warning in melt(genera.agg.16s, id.vars = "Genus", measure.vars =
## c("CSS_Ambient", : The melt generic in data.table has been passed a data.frame
## and will attempt to redirect to the relevant reshape2 method; please note that
## reshape2 is deprecated, and this redirection is now deprecated as well. To
## continue using melt methods from reshape2 while both libraries are attached,
## e.g. melt.list, you can prepend the namespace like reshape2::melt(genera.agg.
## 16s). In the next version, this warning will become an error.
genera.agg.long.16s$Abundance <- ifelse(genera.agg.long.16s$variable == samp.16s[1], genera.agg.long.16s$value / sumsum.16s[1],
                                    ifelse(genera.agg.long.16s$variable == samp.16s[2], genera.agg.long.16s$value / sumsum.16s[2],
                                           ifelse(genera.agg.long.16s$variable == samp.16s[3], genera.agg.long.16s$value / sumsum.16s[3],
                                                  ifelse(genera.agg.long.16s$variable == samp.16s[4], genera.agg.long.16s$value / sumsum.16s[4],"NA" ))))

#call everything below 0.01 (1%) "other_genera"
genera.agg.long.16s$Genus[genera.agg.long.16s$Abundance < 0.01] <- genera.agg.long.16s$Genus["other_genera"] #turn all genera that are below 1% first into "NA" then into "other_genera"
genera.agg.long.16s$Genus <- as.character(genera.agg.long.16s$Genus)
genera.agg.long.16s$Genus[is.na(genera.agg.long.16s$Genus)] <- "other_genera"
genera.agg.long.16s$Abundance <- as.numeric(genera.agg.long.16s$Abundance)

# Checking to make sure it adds to 1 precisely
for (i in 1:length(samp.16s)) {
  print(paste0("For precipitation treatment ", samp.16s[i], " the relative abundance adds to: ", sum(genera.agg.long.16s[genera.agg.long.16s$variable == samp.16s[i], ]$Abundance), "."))
}
## [1] "For precipitation treatment CSS_Ambient the relative abundance adds to: NA."
## [1] "For precipitation treatment CSS_Reduced the relative abundance adds to: NA."
## [1] "For precipitation treatment Grass_Ambient the relative abundance adds to: NA."
## [1] "For precipitation treatment Grass_Reduced the relative abundance adds to: NA."
#create stacked barplot 
genera.agg.long.16s$Genus <- factor(genera.agg.long.16s$Genus)
genera.agg.long.16s$Genus = with(genera.agg.long.16s, reorder(Genus, -Abundance))
#genera.agg.long.16s$Genus <- relevel(genera.agg.long.16s$Genus, 'other_genera') #if want the "other_genera" at the top of the plot.

c36 <- c("#c8c9c9",
"#bda69d",
"#284d50",
"#c1a4a6",
"#344b4a",
"#b2a6bd",
"#58423a",
"#8eb2a9",
"#594152",
"#83adad",
"#56414b",
"#b2aa97",
"#484555",
"#a9ac9f",
"#3a4951",
"#b4a8ab",
"#414849",
"#a1adac",
"#544341",
"#749ca7",
"#5f5740",
"#a98d9e",
"#4b5d49",
"#b4908d",
"#3a586a",
"#896d5f",
"#4e5269",
"#7d8168",
"#665b72",
"#718c7a",
"#6d505d",
"#517368",
"#86656c",
"#46474a",
"#7c819a",
"#967471")

barplot.genus.16s <- ggplot(genera.agg.long.16s, aes(x = fct_relevel(variable, "Grass_Ambient","Grass_Reduced","CSS_Ambient","CSS_Reduced"), y = Abundance, fill = Genus)) +
                     theme_test() +  
                     geom_bar(stat = "identity", width = 0.7) +
                     theme(axis.line.x=element_line(colour = "black"), axis.line.y=element_line(colour = "black")) + 
                     scale_fill_manual(values = c36) +
                     scale_y_continuous(expand = c(0, 0)) +
                     labs(fill = "Genus") +
                     theme(legend.position = "none") + #Change to "right" to display legend at the right of the plot
                     guides(size = "none") +
                     labs(x = "Ecosystem by Water Input", y = "Proportion of Total Bacterial Community")

plot(barplot.genus.16s)
## Warning: Removed 568 rows containing missing values (position_stack).

###Fungi###
taxonomyITS<-read_qza("taxonomyITS.qza", "./tmp", rm=T)
taxonomyITS$uuid
## [1] "e4992fa1-ce2d-4474-b859-f5103b8c08e2"
taxtableITS<-taxonomyITS$data %>% as.tibble() %>% separate(Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species")) #convert the table into a tabular split version
## Warning: Expected 7 pieces. Missing pieces filled with `NA` in 2197 rows [1, 3,
## 4, 5, 6, 7, 8, 10, 11, 13, 16, 17, 20, 21, 22, 23, 24, 25, 26, 27, ...].
taxtableITS
## # A tibble: 4,601 x 9
##    Feature.ID   Kingdom  Phylum  Class  Order Family  Genus  Species  Confidence
##    <fct>        <chr>    <chr>   <chr>  <chr> <chr>   <chr>  <chr>         <dbl>
##  1 7e80c7e5f58… k__Fungi p__Asc… <NA>   <NA>  <NA>    <NA>   <NA>          0.812
##  2 1094975b0bd… k__Fungi p__Asc… c__Le… o__H… f__Hel… g__Ar… s__Arti…      0.952
##  3 3561d8e1d04… k__Fungi p__Asc… c__So… o__G… f__Glo… g__Co… <NA>          0.982
##  4 9995cbaaec4… k__Fungi p__Asc… <NA>   <NA>  <NA>    <NA>   <NA>          0.820
##  5 592482291be… k__Fungi p__Asc… <NA>   <NA>  <NA>    <NA>   <NA>          0.836
##  6 7d8e2be52cd… k__Fungi p__Asc… c__Do… o__C… <NA>    <NA>   <NA>          0.982
##  7 0ca5bba7ce1… k__Fungi p__Asc… c__Do… <NA>  <NA>    <NA>   <NA>          0.836
##  8 53f1a5397d6… k__Plan… p__Chl… c__Tr… o__T… f__Tre… g__Tr… <NA>          0.999
##  9 210acd46015… k__Fungi p__Asc… c__Do… o__P… f__Pha… g__Ch… s__unid…      0.995
## 10 fb069049f2a… k__Fungi <NA>    <NA>   <NA>  <NA>    <NA>   <NA>          0.984
## # … with 4,591 more rows
taxtableITS <- as.data.frame(taxtableITS)

#remove p__ as a precursor to phylum name 
taxtableITS$Phylum <- gsub("p__", "", taxtableITS$Phylum)
taxtableITS$Class <- gsub("c__", "", taxtableITS$Class)
taxtableITS$Order <- gsub("o__", "", taxtableITS$Order)
taxtableITS$Family <- gsub("f__", "", taxtableITS$Family)
taxtableITS$Genus <- gsub("g__", "", taxtableITS$Genus)
taxtableITS$Species <- gsub("s__", "", taxtableITS$Species)

#this will turn all empty cell into NA
taxtableITS <- taxtableITS %>%
  mutate_at(vars(colnames(.)),
            .funs = funs(ifelse(.=="", NA, as.character(.))))

#taxtable was imported using qiime2R, with already separated levels
row.names(taxtableITS) <- taxtableITS$Feature.ID #turn feature IDs into row.names
taxtableITS$Feature.ID <- NULL #remove column "Feature ID"

#merge OTU table with taxonomy table
otu.raref.taxonomy.ITS <- as.data.frame(merge(taxtableITS, rared.otu.ITS.t, by.x = "row.names", by.y = "row.names"))
row.names(otu.raref.taxonomy.ITS) <- otu.raref.taxonomy.ITS$Row.names #move column "Row.names" (first column) to row.names this was created during the merge
otu.raref.taxonomy.ITS$Row.names <- NULL #remove column "Row.names"
otu.raref.taxonomy.ITS$Confidence <- NULL # remove column "Confidence"
otu.raref.taxonomy.ITS.ss <- subset(otu.raref.taxonomy.ITS, select=c(8:195,1:7)) # move Taxa information to the end of all columns
write.table(otu.raref.taxonomy.ITS.ss, "otu_raref_taxonomy.ITS.csv", quote=FALSE, sep=',') #use this table to calculate top number of ESVs  and their sequence contributions.
colSums(otu.raref.taxonomy.ITS.ss[, c(1:188)]) # check if samples are RAREFIED!!
##  14LRNDecember2016ITS  14LRNDecember2017ITS      14LRNJune2017ITS 
##                  1066                  1063                  1066 
##     14LRNMarch2017ITS     14LRNMarch2018ITS 14LRNSeptember2017ITS 
##                  1065                  1064                  1068 
##    14RRXAugust2016ITS  14RRXDecember2016ITS  14RRXDecember2017ITS 
##                  1067                  1060                  1068 
##      14RRXJune2017ITS     14RRXMarch2017ITS     14RRXMarch2018ITS 
##                  1063                  1066                  1067 
## 14RRXSeptember2017ITS    18LXNAugust2016ITS  18LXNDecember2016ITS 
##                  1063                  1065                  1062 
##  18LXNDecember2017ITS      18LXNJune2017ITS     18LXNMarch2017ITS 
##                  1065                  1061                  1068 
##     18LXNMarch2018ITS 18LXNSeptember2017ITS    18RXXAugust2016ITS 
##                  1064                  1063                  1063 
##  18RXXDecember2016ITS  18RXXDecember2017ITS      18RXXJune2017ITS 
##                  1064                  1068                  1064 
##     18RXXMarch2017ITS     18RXXMarch2018ITS 18RXXSeptember2017ITS 
##                  1067                  1058                  1064 
##    20LRXAugust2016ITS  20LRXDecember2016ITS  20LRXDecember2017ITS 
##                  1065                  1063                  1062 
##      20LRXJune2017ITS     20LRXMarch2017ITS     20LRXMarch2018ITS 
##                  1063                  1064                  1063 
## 20LRXSeptember2017ITS  20RRNDecember2016ITS  20RRNDecember2017ITS 
##                  1065                  1068                  1066 
##     20RRNMarch2017ITS     20RRNMarch2018ITS 20RRNSeptember2017ITS 
##                  1067                  1067                  1066 
##    22LXXAugust2016ITS  22LXXDecember2016ITS  22LXXDecember2017ITS 
##                  1068                  1068                  1064 
##      22LXXJune2017ITS     22LXXMarch2017ITS     22LXXMarch2018ITS 
##                  1063                  1066                  1063 
## 22LXXSeptember2017ITS    22RXNAugust2016ITS  22RXNDecember2016ITS 
##                  1064                  1061                  1061 
##  22RXNDecember2017ITS      22RXNJune2017ITS     22RXNMarch2017ITS 
##                  1062                  1066                  1063 
##     22RXNMarch2018ITS 22RXNSeptember2017ITS    25LRXAugust2016ITS 
##                  1063                  1064                  1066 
##  25LRXDecember2016ITS  25LRXDecember2017ITS     25LRXMarch2017ITS 
##                  1061                  1067                  1064 
##     25LRXMarch2018ITS 25LRXSeptember2017ITS    25RRNAugust2016ITS 
##                  1060                  1061                  1066 
##  25RRNDecember2016ITS  25RRNDecember2017ITS      25RRNJune2017ITS 
##                  1061                  1063                  1065 
##     25RRNMarch2017ITS     25RRNMarch2018ITS 25RRNSeptember2017ITS 
##                  1059                  1066                  1065 
##    27LXNAugust2016ITS  27LXNDecember2016ITS  27LXNDecember2017ITS 
##                  1067                  1064                  1065 
##     27LXNMarch2017ITS     27LXNMarch2018ITS 27LXNSeptember2017ITS 
##                  1066                  1063                  1062 
##  27RXXDecember2016ITS      27RXXJune2017ITS     27RXXMarch2017ITS 
##                  1062                  1062                  1064 
##     27RXXMarch2018ITS    32LXNAugust2016ITS  32LXNDecember2016ITS 
##                  1065                  1064                  1065 
##     32LXNMarch2017ITS     32LXNMarch2018ITS 32LXNSeptember2017ITS 
##                  1061                  1062                  1061 
##    32RXXAugust2016ITS  32RXXDecember2016ITS  32RXXDecember2017ITS 
##                  1067                  1063                  1064 
##      32RXXJune2017ITS     32RXXMarch2017ITS     32RXXMarch2018ITS 
##                  1068                  1059                  1066 
## 32RXXSeptember2017ITS    35LRXAugust2016ITS  35LRXDecember2016ITS 
##                  1058                  1064                  1061 
##  35LRXDecember2017ITS      35LRXJune2017ITS     35LRXMarch2017ITS 
##                  1067                  1067                  1066 
##     35LRXMarch2018ITS 35LRXSeptember2017ITS    35RRNAugust2016ITS 
##                  1065                  1064                  1059 
##  35RRNDecember2016ITS  35RRNDecember2017ITS      35RRNJune2017ITS 
##                  1060                  1066                  1059 
##     35RRNMarch2017ITS     35RRNMarch2018ITS 35RRNSeptember2017ITS 
##                  1061                  1065                  1059 
##  45LRXDecember2016ITS  45LRXDecember2017ITS      45LRXJune2017ITS 
##                  1064                  1066                  1062 
##     45LRXMarch2017ITS     45LRXMarch2018ITS 45LRXSeptember2017ITS 
##                  1062                  1065                  1062 
##    45RRNAugust2016ITS  45RRNDecember2016ITS  45RRNDecember2017ITS 
##                  1067                  1066                  1064 
##      45RRNJune2017ITS     45RRNMarch2017ITS     45RRNMarch2018ITS 
##                  1061                  1067                  1062 
## 45RRNSeptember2017ITS    46LXNAugust2016ITS  46LXNDecember2016ITS 
##                  1067                  1059                  1063 
##  46LXNDecember2017ITS      46LXNJune2017ITS     46LXNMarch2017ITS 
##                  1065                  1069                  1062 
##     46LXNMarch2018ITS 46LXNSeptember2017ITS    46RXXAugust2016ITS 
##                  1067                  1063                  1065 
##  46RXXDecember2016ITS  46RXXDecember2017ITS      46RXXJune2017ITS 
##                  1067                  1065                  1061 
##     46RXXMarch2018ITS 46RXXSeptember2017ITS    47LRNAugust2016ITS 
##                  1063                  1067                  1066 
##  47LRNDecember2016ITS  47LRNDecember2017ITS      47LRNJune2017ITS 
##                  1064                  1064                  1063 
##     47LRNMarch2017ITS     47LRNMarch2018ITS 47LRNSeptember2017ITS 
##                  1064                  1068                  1065 
##    47RRXAugust2016ITS  47RRXDecember2016ITS      47RRXJune2017ITS 
##                  1064                  1063                  1061 
##     47RRXMarch2017ITS     47RRXMarch2018ITS 47RRXSeptember2017ITS 
##                  1063                  1066                  1067 
##    48LXXAugust2016ITS  48LXXDecember2016ITS  48LXXDecember2017ITS 
##                  1068                  1064                  1062 
##      48LXXJune2017ITS     48LXXMarch2017ITS     48LXXMarch2018ITS 
##                  1067                  1066                  1072 
##    48RXNAugust2016ITS  48RXNDecember2016ITS  48RXNDecember2017ITS 
##                  1066                  1065                  1068 
##      48RXNJune2017ITS     48RXNMarch2017ITS     48RXNMarch2018ITS 
##                  1065                  1058                  1071 
## 48RXNSeptember2017ITS     4LXXAugust2016ITS   4LXXDecember2016ITS 
##                  1067                  1064                  1066 
##   4LXXDecember2017ITS       4LXXJune2017ITS      4LXXMarch2017ITS 
##                  1064                  1061                  1066 
##  4LXXSeptember2017ITS     4RXNAugust2016ITS   4RXNDecember2017ITS 
##                  1059                  1067                  1064 
##       4RXNJune2017ITS      4RXNMarch2017ITS      4RXNMarch2018ITS 
##                  1063                  1063                  1065 
##  4RXNSeptember2017ITS   5LRNDecember2016ITS       5LRNJune2017ITS 
##                  1063                  1067                  1062 
##      5LRNMarch2017ITS      5LRNMarch2018ITS  5LRNSeptember2017ITS 
##                  1059                  1064                  1061 
##     5RRXAugust2016ITS   5RRXDecember2016ITS   5RRXDecember2017ITS 
##                  1057                  1065                  1064 
##       5RRXJune2017ITS      5RRXMarch2017ITS      5RRXMarch2018ITS 
##                  1071                  1064                  1066 
##  5RRXSeptember2017ITS     7LRXAugust2016ITS   7LRXDecember2016ITS 
##                  1066                  1061                  1065 
##   7LRXDecember2017ITS       7LRXJune2017ITS      7LRXMarch2017ITS 
##                  1066                  1067                  1064 
##      7LRXMarch2018ITS  7LRXSeptember2017ITS   7RRNDecember2016ITS 
##                  1067                  1063                  1066 
##   7RRNDecember2017ITS       7RRNJune2017ITS 
##                  1060                  1067
# Check the number of unique OTUs at different levels
phylum.ITS <- as.vector(unique(otu.raref.taxonomy.ITS.ss$Phylum[!is.na(otu.raref.taxonomy.ITS.ss$Phylum)]))
length(phylum.ITS) #4 phyla
## [1] 4
class.ITS <- as.vector(unique(otu.raref.taxonomy.ITS.ss$Class[!is.na(otu.raref.taxonomy.ITS.ss$Class)]))
length(class.ITS) #21 unique bacterial classes (NA not counted) 
## [1] 21
order.ITS <- as.vector(unique(otu.raref.taxonomy.ITS.ss$Order[!is.na(otu.raref.taxonomy.ITS.ss$Order)]))
length(order.ITS) #67 unique bacterial orders (na not counted) 
## [1] 67
family.ITS <- as.vector(unique(otu.raref.taxonomy.ITS.ss$Family[!is.na(otu.raref.taxonomy.ITS.ss$Family)]))
length(family.ITS) #155 unique families (na not counted)
## [1] 155
genus.ITS <- as.vector(unique(otu.raref.taxonomy.ITS.ss$Genus[!is.na(otu.raref.taxonomy.ITS.ss$Genus)]))
length(genus.ITS) #275 unique bacterial genus' (na not counted) 
## [1] 274
#aggregate taxonomical IDs on a certain level by sample (here at Order level)
genus.ITS.1 <- aggregate(. ~Genus, data = otu.raref.taxonomy.ITS.ss[c(1:(length(names(otu.raref.taxonomy.ITS.ss)) - 7), match("Genus", names(otu.raref.taxonomy.ITS.ss)))], FUN = sum, na.omit=TRUE)

#Order column is turned into rownames and the empty cell is called "unidentified" (o__ gsub into "")
rownames(genus.ITS.1) <- genus.ITS.1$Genus
genus.ITS.1$Genus <- NULL

#turn table into dataframe and transpose
genus.ITS.1 <- as.data.frame(t(genus.ITS.1))
write.xlsx(genus.ITS.1, file = "genera.ITS.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names=TRUE)

# Merge a metadata column here (I just made my own metatdata column)
matching.metadata.gen.ITS <- metadata.ITS.1[c(match(sort(row.names(genus.ITS.1)), sort(row.names(metadata.ITS.1)))),]
merged.meta.gen.ITS <- merge(genus.ITS.1, matching.metadata.gen.ITS, by.x = "row.names", by.y = "row.names")
colnames(merged.meta.gen.ITS)[colnames(merged.meta.gen.ITS)=="Row.names"] <- "SampleID"
merged.meta.gen.ITS <- as.data.frame(merged.meta.gen.ITS)
rownames(merged.meta.gen.ITS) <- merged.meta.gen.ITS$SampleID
merged.meta.gen.ITS$SampleID <- NULL

#aggregate only the X number of the level of taxonomy chosen (from the sequence summaries above) with the metadata selected in this case plot Treatment
genera.agg.ITS <- aggregate(merged.meta.gen.ITS[, 1:275], list(merged.meta.gen.ITS$eco.prec), mean) #genus
## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA

## Warning in mean.default(X[[i]], ...): argument is not numeric or logical:
## returning NA
genera.agg.ITS <- as.data.frame(t(genera.agg.ITS)) # transpose table to have variables as columns
names(genera.agg.ITS) <- as.character(unlist(genera.agg.ITS["Group.1", ])) # transposing causes the treatments named as "Group.1", remove this
genera.agg.ITS <- genera.agg.ITS[-1,] #final removal of "Group.1" 
for(i in c(1:ncol(genera.agg.ITS))) {
  genera.agg.ITS[,i] <- as.numeric(as.character(genera.agg.ITS[,i]))
}
write.xlsx(genera.agg.ITS, file = "plot.genera.ITS.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names=TRUE)

#calculate the number of sequences per level (change the divisor according to above calculation for unique OTUs per level)
totalseqs.ITS <- colSums(genera.agg.ITS) 
avgnumseqs.ITS <- sum(totalseqs.ITS)/275 #142 is the total number of unique genera with NA removed

# CHECK FOR THE MOST ABUNDANT of classifaction level selected HERE!!!!
genera.agg.ITS$Genus <- row.names(genera.agg.ITS)
genera.agg.ITS$Genus <- factor(genera.agg.ITS$Genus)
sumsum.ITS <- colSums(genera.agg.ITS[1:(length(names(genera.agg.ITS))-1)])
samp.ITS <- names(genera.agg.ITS)[1:(length(names(genera.agg.ITS))-1)]

genera.agg.long.ITS <- melt(genera.agg.ITS, id.vars = "Genus", measure.vars = c("CSS_Ambient","CSS_Reduced", "Grass_Ambient","Grass_Reduced"))
## Warning in melt(genera.agg.ITS, id.vars = "Genus", measure.vars =
## c("CSS_Ambient", : The melt generic in data.table has been passed a
## data.frame and will attempt to redirect to the relevant reshape2 method;
## please note that reshape2 is deprecated, and this redirection is now
## deprecated as well. To continue using melt methods from reshape2 while both
## libraries are attached, e.g. melt.list, you can prepend the namespace like
## reshape2::melt(genera.agg.ITS). In the next version, this warning will become an
## error.
genera.agg.long.ITS$Abundance <- ifelse(genera.agg.long.ITS$variable == samp.ITS[1], genera.agg.long.ITS$value / sumsum.ITS[1],
                                    ifelse(genera.agg.long.ITS$variable == samp.ITS[2], genera.agg.long.ITS$value / sumsum.ITS[2],
                                           ifelse(genera.agg.long.ITS$variable == samp.ITS[3], genera.agg.long.ITS$value / sumsum.ITS[3],
                                                  ifelse(genera.agg.long.ITS$variable == samp.ITS[4], genera.agg.long.ITS$value / sumsum.ITS[4],"NA" ))))

#call everything below 0.01 (1%) "other_genera"
genera.agg.long.ITS$Genus[genera.agg.long.ITS$Abundance < 0.01] <- genera.agg.long.ITS$Genus["other_genera"] #turn all genera that are below 1% first into "NA" then into "other_genera"
genera.agg.long.ITS$Genus[ genera.agg.long.ITS$Genus == "unidentified" ] <- NA
genera.agg.long.ITS$Genus <- as.character(genera.agg.long.ITS$Genus)
genera.agg.long.ITS$Genus[is.na(genera.agg.long.ITS$Genus)] <- "other_genera"

genera.agg.long.ITS$Abundance <- as.numeric(genera.agg.long.ITS$Abundance)

# Checking to make sure it adds to 1 precisely
for (i in 1:length(samp.ITS)) {
  print(paste0("For precipitation treatment ", samp.ITS[i], " the relative abundance adds to: ", sum(genera.agg.long.ITS[genera.agg.long.ITS$variable == samp.ITS[i], ]$Abundance), "."))
}
## [1] "For precipitation treatment CSS_Ambient the relative abundance adds to: NA."
## [1] "For precipitation treatment CSS_Reduced the relative abundance adds to: NA."
## [1] "For precipitation treatment Grass_Ambient the relative abundance adds to: NA."
## [1] "For precipitation treatment Grass_Reduced the relative abundance adds to: NA."
#create stacked barplot 
genera.agg.long.ITS$Genus <- factor(genera.agg.long.ITS$Genus)
genera.agg.long.ITS$Genus = with(genera.agg.long.ITS, reorder(Genus, -Abundance))
#genera.agg.long.ITS$Genus <- relevel(genera.agg.long.ITS$Genus, "other_genera") #if want the "other_genera" at the top of the plot.

c18 <- c("#c8c9c9",
         "#264d4c",
         "#bfa4ab",
         "#88abbc",
         "#645b44",
         "#ac9fb8",
         "#736f56",
         "#9692ac",
         "#7c6054",
         "#5f8791",
         "#b7988b",
         "#604b5e",
         "#738b77",
         "#93737e",
         "#5a7c71",
         "#8e8c72",
         "#587789",
         "#676666") 

barplot.genus.ITS <- ggplot(genera.agg.long.ITS, aes(x = fct_relevel(variable, "Grass_Ambient","Grass_Reduced","CSS_Ambient","CSS_Reduced"), y = Abundance, fill = Genus)) +
                     theme_test() +  
                     geom_bar(stat = "identity", width = 0.7) +
                     theme(axis.line.x=element_line(colour = "black"), axis.line.y=element_line(colour = "black")) + 
                     scale_fill_manual(values = c18) +
                     scale_y_continuous(expand = c(0, 0)) +
                     labs(fill = "Genus") +
                     theme(legend.position = "none") + #Change to "right" to display legend at the right of the plot
                     guides(size = "none") +
                     labs(x = "Ecosystem by Water Input", y = "Proportion of Total Fungal Community")

plot(barplot.genus.ITS)
## Warning: Removed 1100 rows containing missing values (position_stack).

###Plant###
plant.func.2015 <- read_excel("plant_barplot.2015.eco.pre.xlsx")
plant.func.2015 <- as.data.frame(plant.func.2015, col.names =   TRUE,   row.names   =   TRUE)
plant.func.2015[1:3 , 1:3]
##   functional_group     community percentage
## 1           g_Bare   css_ambient   2.988300
## 2           g_Bare   css_reduced  25.266378
## 3           g_Bare grass_ambient   2.143758
c07 <- c("#344c60",
         "#baa6ab",
         "#3f4944",
         "#7b9a8b",
         "#85778e",
         "#5c5f47",
         "#82655a")

barplot.plant.2015 <- ggplot(plant.func.2015, aes(x = fct_relevel(community,"grass_ambient","grass_reduced","css_ambient","css_reduced"), y = percentage, fill = functional_group)) +
                                geom_col(aes(fill = factor(functional_group, levels = c("NativeGrassCover", "NonNativeGrassCover","NativeShrubCover","NativeForbCover", "NonNativeForbCover", "Litter",     
                                "Bare")))) +
                                theme_test() +  
                                geom_bar(stat = "identity", width = 0.7) +
                                theme(axis.line.x=element_line(colour = "black"), axis.line.y=element_line(colour = "black")) + 
                                scale_y_continuous(expand = c(0, 0)) +
                                scale_fill_manual(values = c07) +
                                labs(fill = "factor") +
                                theme(legend.position = "none") + #Change to "right" to display legend at the right of the plot
                                guides(size = "right") +
                                labs(x = "Plant Functional Groups by Ecosystem", y = "Proportion Fractional Cover") 

plot(barplot.plant.2015)

###Plot Bacterial, Fungal, and Plant in one figure###
combined.stacked.taxa.genus <- ggarrange(barplot.genus.16s, barplot.genus.ITS, barplot.plant.2015, barplot.plant.2015, labels = c("A.", "B.","C.","D."), ncol = 2, nrow = 2, widths = c(5,5,5,5))
## Warning: Removed 568 rows containing missing values (position_stack).
## Warning: Removed 1100 rows containing missing values (position_stack).
plot(combined.stacked.taxa.genus)

###Richness and Evenness###

###Bacteria###
attach(metadata.16s.1)
Nitrogen<-factor(Nitrogen)
Precipitation <- as.factor(Precipitation)
Block <- as.factor(Block)
Plot <- as.factor(Plot)
Collection_date <- as.factor(Collection_date)
eco.prec <- as.factor(eco.prec)
eco.nit <- as.factor(eco.nit)

bact.aov.richness <- lmer(richness ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date  + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen + ( 1 | Block:Ecosystem), REML = FALSE, data = metadata.16s.1)
## boundary (singular) fit: see ?isSingular
anova(bact.aov.richness, type = "III")
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Block                          2279.8  2279.8     1   192  7.7321 0.0059650 ** 
## Ecosystem                      3545.2  3545.2     1   192 12.0241 0.0006481 ***
## Precipitation                    74.4    74.4     1   192  0.2523 0.6160429    
## Nitrogen                         78.8    78.8     1   192  0.2674 0.6056878    
## Collection_date               12177.5  2029.6     6   192  6.8837 1.234e-06 ***
## Ecosystem:Precipitation         526.7   526.7     1   192  1.7863 0.1829637    
## Ecosystem:Nitrogen               91.9    91.9     1   192  0.3117 0.5772818    
## Ecosystem:Collection_date      2552.8   425.5     6   192  1.4430 0.2001184    
## Precipitation:Collection_date  3396.3   566.1     6   192  1.9199 0.0794684 .  
## Nitrogen:Collection_date       1137.2   189.5     6   192  0.6428 0.6958267    
## Precipitation:Nitrogen         1373.4  1373.4     1   192  4.6580 0.0321501 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bact.eco.prec.rich <- ggplot(data = bact.aov.richness, aes(x = factor(eco.prec, levels = c("Grass_Ambient","Grass_Reduced","CSS_Ambient", "CSS_Reduced")), y = metadata.16s.1$richness, fill = metadata.16s.1$Precipitation)) + 
                      geom_boxplot() +
                      geom_dotplot(binaxis='y', stackdir='center',
                      position=position_dodge(0.75), dotsize = 0.65) +
                      scale_fill_manual(values=c('#408080', '#FF8040')) +
                      labs(title = 'Bacterial Richness by Ecosystem', x = 'Ecosystem', y = 'Richness', fill = 'Precipitation') +
                      theme_classic() +
                      theme(panel.grid.major = element_blank(), 
                      panel.grid.minor = element_blank(), 
                      panel.background = element_blank(),
                      panel.border = element_blank())

plot(bact.eco.prec.rich)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

bact.eco.prec.tuk.rich <- TukeyHSD(aov(formula = metadata.16s.1$richness ~ as.character(metadata.16s.1$eco.prec)))
bact.eco.prec.tuk.rich
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = metadata.16s.1$richness ~ as.character(metadata.16s.1$eco.prec))
## 
## $`as.character(metadata.16s.1$eco.prec)`
##                                   diff        lwr       upr     p adj
## CSS_Reduced-CSS_Ambient     -4.9400000 -15.657512  5.777512 0.6308928
## Grass_Ambient-CSS_Ambient   -7.7208333 -18.549410  3.107744 0.2542266
## Grass_Reduced-CSS_Ambient   -4.8363636 -15.913221  6.240494 0.6703734
## Grass_Ambient-CSS_Reduced   -2.7808333 -13.609410  8.047744 0.9098615
## Grass_Reduced-CSS_Reduced    0.1036364 -10.973221 11.180494 0.9999949
## Grass_Reduced-Grass_Ambient  2.8844697  -8.299885 14.068824 0.9088281
bact.eco.nit.rich <- ggplot(data = bact.aov.richness, aes(x = factor(eco.nit, levels = c("Grass_Ambient","Grass_Added","CSS_Ambient", "CSS_Added")), y = metadata.16s.1$richness, fill = metadata.16s.1$Nitrogen)) + 
                     geom_boxplot() +
                     geom_dotplot(binaxis='y', stackdir='center',
                     position=position_dodge(0.75), dotsize = 0.65) +
                     scale_fill_manual(values=c('#b7a28c','#408080')) +
                     labs(title = 'Bacterial Richness by Ecosystem', x = 'Ecosystem', y = 'Richness', fill = 'Nitrogen') +
                     theme_classic() +
                     theme(panel.grid.major = element_blank(), 
                     panel.grid.minor = element_blank(), 
                     panel.background = element_blank(),
                     panel.border = element_blank())

plot(bact.eco.nit.rich)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

bact.aov.shannon <- lmer(shannon ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date  + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen + ( 1 | Block:Ecosystem), REML = FALSE, data = metadata.16s.1)
anova(bact.aov.shannon, type = "III")
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)
## Block                         0.8067 0.80671     1   7.682  3.5760  0.096799
## Ecosystem                     1.9936 1.99360     1   7.609  8.8373  0.018832
## Precipitation                 0.2835 0.28350     1 184.221  1.2567  0.263732
## Nitrogen                      0.1202 0.12023     1 186.131  0.5330  0.466284
## Collection_date               9.5941 1.59902     6 184.997  7.0882 8.248e-07
## Ecosystem:Precipitation       0.0268 0.02678     1 184.032  0.1187  0.730827
## Ecosystem:Nitrogen            0.1979 0.19785     1 185.400  0.8771  0.350228
## Ecosystem:Collection_date     2.2662 0.37771     6 185.186  1.6743  0.129466
## Precipitation:Collection_date 2.9398 0.48997     6 184.617  2.1719  0.047575
## Nitrogen:Collection_date      1.2565 0.20941     6 185.166  0.9283  0.475847
## Precipitation:Nitrogen        1.8126 1.81260     1 184.929  8.0349  0.005099
##                                  
## Block                         .  
## Ecosystem                     *  
## Precipitation                    
## Nitrogen                         
## Collection_date               ***
## Ecosystem:Precipitation          
## Ecosystem:Nitrogen               
## Ecosystem:Collection_date        
## Precipitation:Collection_date *  
## Nitrogen:Collection_date         
## Precipitation:Nitrogen        ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
bact.eco.prec.shannon <- ggplot(data = bact.aov.shannon, aes(x = factor(eco.prec, levels = c("Grass_Ambient","Grass_Reduced","CSS_Ambient", "CSS_Reduced")), y = metadata.16s.1$shannon, fill = metadata.16s.1$Precipitation)) + 
                         geom_boxplot() +
                         geom_dotplot(binaxis='y', stackdir='center',
                         position=position_dodge(0.75), dotsize = 0.65) +
                         scale_fill_manual(values=c('#408080', '#FF8040')) +
                         labs(title = 'Bacterial Evenness by Ecosystem', x = 'Ecosystem', y = 'Evenness', fill = 'Precipitation') +
                         theme_classic() +
                         #theme(legend.position = "none") +
                         theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())

bact.eco.prec.tuk.shannon <- TukeyHSD(aov(formula = metadata.16s.1$shannon ~ as.character(metadata.16s.1$eco.prec)))
bact.eco.prec.tuk.shannon
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = metadata.16s.1$shannon ~ as.character(metadata.16s.1$eco.prec))
## 
## $`as.character(metadata.16s.1$eco.prec)`
##                                    diff        lwr        upr     p adj
## CSS_Reduced-CSS_Ambient     -0.06079776 -0.3615575 0.23996197 0.9532127
## Grass_Ambient-CSS_Ambient   -0.16421605 -0.4680925 0.13966044 0.5004130
## Grass_Reduced-CSS_Ambient   -0.23491663 -0.5457605 0.07592721 0.2073915
## Grass_Ambient-CSS_Reduced   -0.10341829 -0.4072948 0.20045820 0.8141064
## Grass_Reduced-CSS_Reduced   -0.17411887 -0.4849627 0.13672497 0.4687036
## Grass_Reduced-Grass_Ambient -0.07070058 -0.3845611 0.24315991 0.9368327
bact.eco.nit.shannon <- ggplot(data = bact.aov.shannon, aes(x = factor(eco.nit, levels = c("Grass_Ambient","Grass_Added","CSS_Ambient", "CSS_Added")), y = metadata.16s.1$shannon, fill = metadata.16s.1$Nitrogen)) + 
                        geom_boxplot() +
                        geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.75), dotsize = 0.65) +
                        scale_fill_manual(values=c('#b7a28c','#408080')) +
                        labs(title = 'Bacterial Evenness by Ecosystem', x = 'Ecosystem', y = 'Evenness', fill = 'Nitrogen') +
                        theme_classic() +
                        #theme(legend.position = "none") +
                        theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())

plot(bact.eco.nit.shannon)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

detach(metadata.16s.1)

###Fungi###
attach(metadata.ITS.1)
## The following objects are masked _by_ .GlobalEnv:
## 
##     Block, Collection_date, eco.nit, eco.prec, Nitrogen, Plot,
##     Precipitation
Nitrogen<-factor(Nitrogen)
Precipitation <- as.factor(Precipitation)
Block <- as.factor(Block)
Plot <- as.factor(Plot)
Collection_date <- as.factor(Collection_date)
eco.prec <- as.factor(eco.prec)
eco.nit <- as.factor(eco.nit)

fung.aov.richness <- lmer(richness ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date  + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen + ( 1 | Block:Ecosystem), REML = FALSE, data = metadata.ITS.1)
## boundary (singular) fit: see ?isSingular
anova(fung.aov.richness, type = "III")
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Block                            59.9    59.9     1   203  0.0634 0.801500   
## Ecosystem                      8863.5  8863.5     1   203  9.3822 0.002488 **
## Precipitation                  3175.5  3175.5     1   203  3.3613 0.068209 . 
## Nitrogen                         89.9    89.9     1   203  0.0951 0.758053   
## Collection_date                8574.7  1429.1     6   203  1.5127 0.175484   
## Ecosystem:Precipitation        4453.3  4453.3     1   203  4.7139 0.031080 * 
## Ecosystem:Nitrogen             1304.9  1304.9     1   203  1.3813 0.241258   
## Ecosystem:Collection_date     11942.0  1990.3     6   203  2.1068 0.053957 . 
## Precipitation:Collection_date  7921.5  1320.2     6   203  1.3975 0.217150   
## Nitrogen:Collection_date      12097.0  2016.2     6   203  2.1341 0.050978 . 
## Precipitation:Nitrogen          748.6   748.6     1   203  0.7924 0.374438   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fung.eco.prec.rich <- ggplot(data = fung.aov.richness, aes(x = factor(eco.prec, levels = c("Grass_Ambient","Grass_Reduced","CSS_Ambient", "CSS_Reduced")), y = metadata.ITS.1$richness, fill = metadata.ITS.1$Precipitation)) + 
                      geom_boxplot() +
                      geom_dotplot(binaxis='y', stackdir='center',
                      position=position_dodge(0.75), dotsize = 0.65) +
                      scale_fill_manual(values=c('#408080', '#FF8040')) +
                      labs(title = 'Fungal Richness by Ecosystem', x = 'Ecosystem', y = 'Richness', fill = 'Precipitation') +
                      theme_classic() +
                      theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())

plot(fung.eco.prec.rich)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

fung.eco.prec.tuk.rich <- TukeyHSD(aov(formula = metadata.ITS.1$richness ~ as.character(metadata.ITS.1$eco.prec)))
fung.eco.prec.tuk.rich
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = metadata.ITS.1$richness ~ as.character(metadata.ITS.1$eco.prec))
## 
## $`as.character(metadata.ITS.1$eco.prec)`
##                                   diff       lwr         upr     p adj
## CSS_Reduced-CSS_Ambient     -18.877358 -36.41638  -1.3383342 0.0294106
## Grass_Ambient-CSS_Ambient   -34.726415 -52.26544 -17.1873908 0.0000041
## Grass_Reduced-CSS_Ambient   -36.071429 -53.94744 -18.1954220 0.0000026
## Grass_Ambient-CSS_Reduced   -15.849057 -32.94844   1.2503233 0.0800862
## Grass_Reduced-CSS_Reduced   -17.194070 -34.63893   0.2507867 0.0550132
## Grass_Reduced-Grass_Ambient  -1.345013 -18.78987  16.0998433 0.9971695
fung.eco.nit.rich <- ggplot(data = fung.aov.richness, aes(x = factor(eco.nit, levels = c("Grass_Ambient","Grass_Added","CSS_Ambient", "CSS_Added")), y = metadata.ITS.1$richness, fill = metadata.ITS.1$Nitrogen)) + 
                     geom_boxplot() +
                     geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.75), dotsize = 0.65) +
                     scale_fill_manual(values=c('#b7a28c','#408080')) +
                     labs(title = 'Fungal Richness by Ecosystem', x = 'Ecosystem', y = 'Richness', fill = 'Nitrogen') +
                     theme_classic() +
                     theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())

plot(fung.eco.nit.rich)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

fung.aov.shannon <- lmer(shannon ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date  + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen + ( 1 | Block:Ecosystem), REML = FALSE, data = metadata.ITS.1)
## boundary (singular) fit: see ?isSingular
anova(fung.aov.shannon, type = "III")
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## Block                         0.0026 0.00256     1   203  0.0152 0.902014   
## Ecosystem                     0.9103 0.91032     1   203  5.4137 0.020965 * 
## Precipitation                 0.2116 0.21159     1   203  1.2584 0.263286   
## Nitrogen                      0.0003 0.00033     1   203  0.0020 0.964757   
## Collection_date               2.6120 0.43533     6   203  2.5889 0.019388 * 
## Ecosystem:Precipitation       0.7041 0.70406     1   203  4.1871 0.042021 * 
## Ecosystem:Nitrogen            0.3806 0.38064     1   203  2.2637 0.133992   
## Ecosystem:Collection_date     3.2525 0.54208     6   203  3.2238 0.004786 **
## Precipitation:Collection_date 1.8596 0.30994     6   203  1.8432 0.092353 . 
## Nitrogen:Collection_date      1.9503 0.32506     6   203  1.9331 0.077052 . 
## Precipitation:Nitrogen        0.0168 0.01680     1   203  0.0999 0.752270   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fung.eco.prec.shannon <- ggplot(data = fung.aov.shannon, aes(x = factor(eco.prec, levels = c("Grass_Ambient","Grass_Reduced","CSS_Ambient", "CSS_Reduced")), y = metadata.ITS.1$shannon, fill = metadata.ITS.1$Precipitation)) + 
                         geom_boxplot() +
                         geom_dotplot(binaxis='y', stackdir='center',
                         position=position_dodge(0.75), dotsize = 0.65) +
                         scale_fill_manual(values=c('#408080', '#FF8040')) +
                         labs(title = 'Fungal Evenness by Ecosystem', x = 'Ecosystem', y = 'Evenness', fill = 'Precipitation') +
                         theme_classic() +
                         theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())

plot(fung.eco.prec.shannon)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

fung.eco.prec.tuk.shannon <- TukeyHSD(aov(formula = metadata.ITS.1$shannon ~ as.character(metadata.ITS.1$eco.prec)))
fung.eco.prec.tuk.shannon
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = metadata.ITS.1$shannon ~ as.character(metadata.ITS.1$eco.prec))
## 
## $`as.character(metadata.ITS.1$eco.prec)`
##                                    diff        lwr         upr     p adj
## CSS_Reduced-CSS_Ambient     -0.20575542 -0.4467623  0.03525147 0.1235261
## Grass_Ambient-CSS_Ambient   -0.41592122 -0.6569281 -0.17491433 0.0000766
## Grass_Reduced-CSS_Ambient   -0.39190067 -0.6375381 -0.14626325 0.0003052
## Grass_Ambient-CSS_Reduced   -0.21016580 -0.4451315  0.02479986 0.0975101
## Grass_Reduced-CSS_Reduced   -0.18614525 -0.4258582  0.05356766 0.1870450
## Grass_Reduced-Grass_Ambient  0.02402055 -0.2156924  0.26373346 0.9938505
fung.eco.nit.shannon <- ggplot(data = fung.aov.shannon, aes(x = factor(eco.nit, levels = c("Grass_Ambient","Grass_Added","CSS_Ambient", "CSS_Added")), y = metadata.ITS.1$shannon, fill = metadata.ITS.1$Nitrogen)) + 
                        geom_boxplot() +
                        geom_dotplot(binaxis='y', stackdir='center', position=position_dodge(0.75), dotsize = 0.65) +
                        scale_fill_manual(values=c('#b7a28c','#408080')) +
                        labs(title = 'Fungal Evenness by Ecosystem', x = 'Ecosystem', y = 'Evenness', fill = 'Nitrogen') +
                        theme_classic() +
                        theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), panel.border = element_blank())


plot(fung.eco.nit.shannon)
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.

combined_rich_shan_bact_fun.precipitation <- ggarrange(bact.eco.prec.rich, bact.eco.prec.shannon, fung.eco.prec.rich, fung.eco.prec.shannon, labels = c("A.", "B.", "C.", "D."), ncol = 2, nrow = 2, widths = c(5,5)) # to remove legends and titles add "#" next to labs options of the individual plots.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
plot(combined_rich_shan_bact_fun.precipitation)

combined_rich_shan_bact_fun.nitrogen <- ggarrange(bact.eco.nit.rich, bact.eco.nit.shannon, fung.eco.nit.rich, fung.eco.nit.shannon, labels = c("A.", "B.", "C.", "D."), ncol = 2, nrow = 2, widths = c(5,5))
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bindot()` using `bins = 30`. Pick better value with `binwidth`.
plot(combined_rich_shan_bact_fun.nitrogen)

detach(metadata.ITS.1)
###Bacteria###
median.avg.dist.16s <- avgdist(otu.clean.16s, sample = 1090, iterations = 1000, meanfun = median, dmethod = "bray")
## Warning in avgdist(otu.clean.16s, sample = 1090, iterations = 1000, meanfun
## = median, : The following sampling units were removed because they were below
## sampling depth: 14LRNDecember201716S, 14LRNSeptember201716S, 20RRNJune201716S,
## 20RRNSeptember201716S, 22LXXJune201716S, 22LXXSeptember201716S,
## 25LRXJune201716S, 32LXNDecember201716S, 35RRNJune201716S, 45LRXDecember201716S,
## 46RXXDecember201716S, 47LRNDecember201716S, 47RRXDecember201716S,
## 48LXXSeptember201716S, 48RXNDecember201716S, 4LXXAugust201616S, 4LXXJune201716S,
## 5LRNJune201716S, 8RXNJune201716S
bray.dist.16s <- as.data.frame(as.matrix(median.avg.dist.16s))
write.xlsx2(bray.dist.16s, file = "bray.dist.16s.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names = TRUE, append = FALSE)
NMDS.16s <- metaMDS(bray.dist.16s, distance = "bray", k = 3, noshare = FALSE, trymax = 400, trace = TRUE)
## Run 0 stress 0.1473617 
## Run 1 stress 0.1470132 
## ... New best solution
## ... Procrustes: rmse 0.01099409  max resid 0.109775 
## Run 2 stress 0.1471217 
## ... Procrustes: rmse 0.007527532  max resid 0.08744722 
## Run 3 stress 0.147153 
## ... Procrustes: rmse 0.009782592  max resid 0.1102538 
## Run 4 stress 0.147169 
## ... Procrustes: rmse 0.008878901  max resid 0.08772956 
## Run 5 stress 0.1471381 
## ... Procrustes: rmse 0.008596456  max resid 0.1099102 
## Run 6 stress 0.1474069 
## ... Procrustes: rmse 0.00947325  max resid 0.08680287 
## Run 7 stress 0.1471773 
## ... Procrustes: rmse 0.008628991  max resid 0.08740841 
## Run 8 stress 0.1470301 
## ... Procrustes: rmse 0.002391889  max resid 0.02444052 
## Run 9 stress 0.1470599 
## ... Procrustes: rmse 0.005266206  max resid 0.05515962 
## Run 10 stress 0.1472593 
## ... Procrustes: rmse 0.0115494  max resid 0.1112863 
## Run 11 stress 0.1471472 
## ... Procrustes: rmse 0.008735791  max resid 0.1100061 
## Run 12 stress 0.147335 
## ... Procrustes: rmse 0.01178672  max resid 0.109496 
## Run 13 stress 0.1471293 
## ... Procrustes: rmse 0.007341031  max resid 0.08740262 
## Run 14 stress 0.1470149 
## ... Procrustes: rmse 0.0001713161  max resid 0.001736122 
## ... Similar to previous best
## Run 15 stress 0.1470132 
## ... New best solution
## ... Procrustes: rmse 0.00212278  max resid 0.02365873 
## Run 16 stress 0.147158 
## ... Procrustes: rmse 0.008449116  max resid 0.1116341 
## Run 17 stress 0.1472593 
## ... Procrustes: rmse 0.01153399  max resid 0.1113063 
## Run 18 stress 0.1472745 
## ... Procrustes: rmse 0.01167785  max resid 0.1134128 
## Run 19 stress 0.1472637 
## ... Procrustes: rmse 0.01155497  max resid 0.111502 
## Run 20 stress 0.146991 
## ... New best solution
## ... Procrustes: rmse 0.001827964  max resid 0.0203683 
## Run 21 stress 0.147257 
## ... Procrustes: rmse 0.011816  max resid 0.1111026 
## Run 22 stress 0.14701 
## ... Procrustes: rmse 0.00176931  max resid 0.02012093 
## Run 23 stress 0.147356 
## ... Procrustes: rmse 0.009666349  max resid 0.08922747 
## Run 24 stress 0.1471465 
## ... Procrustes: rmse 0.006709326  max resid 0.0868903 
## Run 25 stress 0.1472624 
## ... Procrustes: rmse 0.01203293  max resid 0.1115047 
## Run 26 stress 0.147162 
## ... Procrustes: rmse 0.01013879  max resid 0.1127369 
## Run 27 stress 0.1470135 
## ... Procrustes: rmse 0.002892814  max resid 0.02483455 
## Run 28 stress 0.1472321 
## ... Procrustes: rmse 0.01145066  max resid 0.111831 
## Run 29 stress 0.1469862 
## ... New best solution
## ... Procrustes: rmse 0.002212889  max resid 0.02459004 
## Run 30 stress 0.1472525 
## ... Procrustes: rmse 0.01189783  max resid 0.1113853 
## Run 31 stress 0.1472532 
## ... Procrustes: rmse 0.01189261  max resid 0.1113517 
## Run 32 stress 0.1471686 
## ... Procrustes: rmse 0.008376918  max resid 0.08749068 
## Run 33 stress 0.1470127 
## ... Procrustes: rmse 0.002762964  max resid 0.02300611 
## Run 34 stress 0.1471424 
## ... Procrustes: rmse 0.009824468  max resid 0.1099773 
## Run 35 stress 0.1470191 
## ... Procrustes: rmse 0.003028171  max resid 0.02300531 
## Run 36 stress 0.14707 
## ... Procrustes: rmse 0.00526875  max resid 0.05637375 
## Run 37 stress 0.1472604 
## ... Procrustes: rmse 0.01187004  max resid 0.1114016 
## Run 38 stress 0.1471355 
## ... Procrustes: rmse 0.006614468  max resid 0.08691199 
## Run 39 stress 0.1471404 
## ... Procrustes: rmse 0.0090955  max resid 0.1102109 
## Run 40 stress 0.1471245 
## ... Procrustes: rmse 0.007358757  max resid 0.08761498 
## Run 41 stress 0.1471608 
## ... Procrustes: rmse 0.008310955  max resid 0.1103223 
## Run 42 stress 0.1473535 
## ... Procrustes: rmse 0.01084271  max resid 0.1099971 
## Run 43 stress 0.1472248 
## ... Procrustes: rmse 0.01112402  max resid 0.1112418 
## Run 44 stress 0.1471302 
## ... Procrustes: rmse 0.007576645  max resid 0.08783906 
## Run 45 stress 0.1470569 
## ... Procrustes: rmse 0.0046187  max resid 0.05439872 
## Run 46 stress 0.1469986 
## ... Procrustes: rmse 0.001149863  max resid 0.01225147 
## Run 47 stress 0.1471684 
## ... Procrustes: rmse 0.008967158  max resid 0.1072389 
## Run 48 stress 0.1472329 
## ... Procrustes: rmse 0.01143397  max resid 0.1111328 
## Run 49 stress 0.1472266 
## ... Procrustes: rmse 0.01131513  max resid 0.1112707 
## Run 50 stress 0.1472576 
## ... Procrustes: rmse 0.01204327  max resid 0.1115062 
## Run 51 stress 0.1470195 
## ... Procrustes: rmse 0.002260496  max resid 0.01956336 
## Run 52 stress 0.1471483 
## ... Procrustes: rmse 0.008977734  max resid 0.1106712 
## Run 53 stress 0.1473286 
## ... Procrustes: rmse 0.006317869  max resid 0.06006143 
## Run 54 stress 0.1471239 
## ... Procrustes: rmse 0.006817903  max resid 0.08728222 
## Run 55 stress 0.147125 
## ... Procrustes: rmse 0.007219689  max resid 0.08716617 
## Run 56 stress 0.1470168 
## ... Procrustes: rmse 0.002097874  max resid 0.02017475 
## Run 57 stress 0.1471868 
## ... Procrustes: rmse 0.00858156  max resid 0.08732222 
## Run 58 stress 0.1472545 
## ... Procrustes: rmse 0.01195251  max resid 0.1113325 
## Run 59 stress 0.1471491 
## ... Procrustes: rmse 0.009831946  max resid 0.1101909 
## Run 60 stress 0.1471254 
## ... Procrustes: rmse 0.006842557  max resid 0.08680285 
## Run 61 stress 0.1472416 
## ... Procrustes: rmse 0.008796869  max resid 0.08692925 
## Run 62 stress 0.1470671 
## ... Procrustes: rmse 0.005111794  max resid 0.05578186 
## Run 63 stress 0.1471321 
## ... Procrustes: rmse 0.00659307  max resid 0.08692297 
## Run 64 stress 0.1473484 
## ... Procrustes: rmse 0.01087259  max resid 0.1100298 
## Run 65 stress 0.1472363 
## ... Procrustes: rmse 0.01149527  max resid 0.1110132 
## Run 66 stress 0.1471272 
## ... Procrustes: rmse 0.007194208  max resid 0.08774139 
## Run 67 stress 0.1472331 
## ... Procrustes: rmse 0.01144492  max resid 0.111155 
## Run 68 stress 0.1472372 
## ... Procrustes: rmse 0.01099134  max resid 0.1114732 
## Run 69 stress 0.1470017 
## ... Procrustes: rmse 0.00258371  max resid 0.02365188 
## Run 70 stress 0.1470111 
## ... Procrustes: rmse 0.001886286  max resid 0.01936161 
## Run 71 stress 0.1472254 
## ... Procrustes: rmse 0.0112731  max resid 0.1112213 
## Run 72 stress 0.1473473 
## ... Procrustes: rmse 0.01099156  max resid 0.109975 
## Run 73 stress 0.1472684 
## ... Procrustes: rmse 0.01199933  max resid 0.1115151 
## Run 74 stress 0.1471422 
## ... Procrustes: rmse 0.006670685  max resid 0.0869103 
## Run 75 stress 0.1473327 
## ... Procrustes: rmse 0.01174141  max resid 0.1098759 
## Run 76 stress 0.1472306 
## ... Procrustes: rmse 0.01083578  max resid 0.1113879 
## Run 77 stress 0.1471704 
## ... Procrustes: rmse 0.008270334  max resid 0.08771603 
## Run 78 stress 0.1473392 
## ... Procrustes: rmse 0.01176678  max resid 0.1095753 
## Run 79 stress 0.1471894 
## ... Procrustes: rmse 0.007127218  max resid 0.09445938 
## Run 80 stress 0.1473505 
## ... Procrustes: rmse 0.01081933  max resid 0.1099293 
## Run 81 stress 0.1470164 
## ... Procrustes: rmse 0.002854988  max resid 0.02567395 
## Run 82 stress 0.1471378 
## ... Procrustes: rmse 0.007406443  max resid 0.0878771 
## Run 83 stress 0.1473474 
## ... Procrustes: rmse 0.0108633  max resid 0.1099446 
## Run 84 stress 0.1472608 
## ... Procrustes: rmse 0.01203316  max resid 0.1111735 
## Run 85 stress 0.1470773 
## ... Procrustes: rmse 0.004053615  max resid 0.04704992 
## Run 86 stress 0.1474857 
## ... Procrustes: rmse 0.008960734  max resid 0.08628311 
## Run 87 stress 0.1470152 
## ... Procrustes: rmse 0.002843881  max resid 0.02242059 
## Run 88 stress 0.1469966 
## ... Procrustes: rmse 0.002276933  max resid 0.02330802 
## Run 89 stress 0.1472297 
## ... Procrustes: rmse 0.01137699  max resid 0.1111843 
## Run 90 stress 0.1471237 
## ... Procrustes: rmse 0.006684047  max resid 0.08705333 
## Run 91 stress 0.1471449 
## ... Procrustes: rmse 0.009775783  max resid 0.1101108 
## Run 92 stress 0.1472231 
## ... Procrustes: rmse 0.01110888  max resid 0.1112124 
## Run 93 stress 0.1472235 
## ... Procrustes: rmse 0.01115106  max resid 0.1112283 
## Run 94 stress 0.1470704 
## ... Procrustes: rmse 0.003671529  max resid 0.04225642 
## Run 95 stress 0.1471443 
## ... Procrustes: rmse 0.01001536  max resid 0.1104978 
## Run 96 stress 0.1473516 
## ... Procrustes: rmse 0.01084375  max resid 0.1100857 
## Run 97 stress 0.1471432 
## ... Procrustes: rmse 0.009895986  max resid 0.1096512 
## Run 98 stress 0.147237 
## ... Procrustes: rmse 0.01144252  max resid 0.111118 
## Run 99 stress 0.1470921 
## ... Procrustes: rmse 0.005504368  max resid 0.05578434 
## Run 100 stress 0.1471374 
## ... Procrustes: rmse 0.006686398  max resid 0.08706045 
## Run 101 stress 0.1471497 
## ... Procrustes: rmse 0.01004118  max resid 0.1105268 
## Run 102 stress 0.1471222 
## ... Procrustes: rmse 0.006576593  max resid 0.08696867 
## Run 103 stress 0.1473412 
## ... Procrustes: rmse 0.01173908  max resid 0.1102461 
## Run 104 stress 0.1471212 
## ... Procrustes: rmse 0.00737984  max resid 0.08762399 
## Run 105 stress 0.1471201 
## ... Procrustes: rmse 0.007071948  max resid 0.08758696 
## Run 106 stress 0.1472333 
## ... Procrustes: rmse 0.01026069  max resid 0.1110523 
## Run 107 stress 0.1473655 
## ... Procrustes: rmse 0.01164181  max resid 0.1051452 
## Run 108 stress 0.1472238 
## ... Procrustes: rmse 0.01117147  max resid 0.1112131 
## Run 109 stress 0.1472509 
## ... Procrustes: rmse 0.01177104  max resid 0.1113409 
## Run 110 stress 0.1471873 
## ... Procrustes: rmse 0.008487483  max resid 0.08722651 
## Run 111 stress 0.1473536 
## ... Procrustes: rmse 0.01094308  max resid 0.1098447 
## Run 112 stress 0.1472272 
## ... Procrustes: rmse 0.01132252  max resid 0.1112118 
## Run 113 stress 0.1469921 
## ... Procrustes: rmse 0.002164316  max resid 0.02400088 
## Run 114 stress 0.147131 
## ... Procrustes: rmse 0.006668046  max resid 0.08717988 
## Run 115 stress 0.1471793 
## ... Procrustes: rmse 0.008121136  max resid 0.08704525 
## Run 116 stress 0.1471438 
## ... Procrustes: rmse 0.009907185  max resid 0.1100165 
## Run 117 stress 0.1473549 
## ... Procrustes: rmse 0.0108413  max resid 0.1099318 
## Run 118 stress 0.1470648 
## ... Procrustes: rmse 0.004396068  max resid 0.05314428 
## Run 119 stress 0.1472324 
## ... Procrustes: rmse 0.01143582  max resid 0.1111719 
## Run 120 stress 0.1472248 
## ... Procrustes: rmse 0.01123424  max resid 0.1112131 
## Run 121 stress 0.1471215 
## ... Procrustes: rmse 0.006601118  max resid 0.08702806 
## Run 122 stress 0.1472351 
## ... Procrustes: rmse 0.01135566  max resid 0.1111677 
## Run 123 stress 0.1472358 
## ... Procrustes: rmse 0.01154947  max resid 0.1111586 
## Run 124 stress 0.14725 
## ... Procrustes: rmse 0.01137833  max resid 0.1111902 
## Run 125 stress 0.1473318 
## ... Procrustes: rmse 0.01176556  max resid 0.1100493 
## Run 126 stress 0.1469894 
## ... Procrustes: rmse 0.0006347364  max resid 0.004565108 
## ... Similar to previous best
## *** Solution reached
stressplot(NMDS.16s)

coordinates.16s <- data.frame(NMDS.16s$points[,1:3])
nmds.metadata.16s <- merge(coordinates.16s, metadata.16s.1, by.x = "row.names", by.y = "row.names")
nmds.metadata.16s <- data.frame(nmds.metadata.16s, col.names    =   TRUE,   row.names   =   TRUE)


bacteria.nmds.plot <-  ggplot(data = nmds.metadata.16s) +
                       aes(x = MDS1, y = MDS3, color = as.factor(Precipitation)) +
                       geom_point(aes(shape = as.factor(Ecosystem), size = as.factor(Ecosystem))) +
                       scale_color_manual(values = c('#408080', '#FF8040')) +
                       scale_shape_manual(values = c(17,16)) +
                       scale_size_manual(values = c(5,5)) +
                       geom_text(label = as.factor(nmds.metadata.16s$eco.prec), size = 1) +
                       theme_minimal() + 
                       theme(legend.position = "none") +
                       labs(col = "Precipitation", shape = "Ecosystem") + 
                       theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                       panel.background = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=1), axis.text.x=element_blank(), 
                       axis.ticks.x=element_blank(),axis.text.y=element_blank(),  
                       axis.ticks.y=element_blank())

plot(bacteria.nmds.plot)

PERMANOVA.16s <- adonis(bray.dist.16s ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen, data = nmds.metadata.16s, strata = nmds.metadata.16s$Block, permutations = 999)

PERMANOVA.16s
## 
## Call:
## adonis(formula = bray.dist.16s ~ Block + Ecosystem * Precipitation +      Ecosystem * Nitrogen + Ecosystem * Collection_date + Collection_date *      Precipitation + Collection_date * Nitrogen + Precipitation *      Nitrogen, data = nmds.metadata.16s, permutations = 999, strata = nmds.metadata.16s$Block) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                                Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Block                           1     3.206  3.2060 11.6324 0.04560  0.001 ***
## Ecosystem                       1     4.019  4.0194 14.5836 0.05717  0.001 ***
## Precipitation                   1     1.795  1.7948  6.5119 0.02553  0.001 ***
## Nitrogen                        1     0.632  0.6316  2.2916 0.00898  0.001 ***
## Collection_date                 6     7.456  1.2427  4.5090 0.10606  0.001 ***
## Ecosystem:Precipitation         1     1.283  1.2833  4.6563 0.01825  0.001 ***
## Ecosystem:Nitrogen              1     0.511  0.5109  1.8537 0.00727  0.003 ** 
## Ecosystem:Collection_date       6     2.805  0.4676  1.6965 0.03990  0.001 ***
## Precipitation:Collection_date   6     2.310  0.3851  1.3971 0.03286  0.001 ***
## Nitrogen:Collection_date        6     1.666  0.2777  1.0075 0.02370  0.337    
## Precipitation:Nitrogen          1     0.524  0.5241  1.9015 0.00745  0.004 ** 
## Residuals                     160    44.098  0.2756         0.62723           
## Total                         191    70.306                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###Fungi###
median.avg.dist.ITS <- avgdist(otu.clean.ITS, sample = 1064, iterations = 1000, meanfun = median, dmethod = "bray")
## Warning in avgdist(otu.clean.ITS, sample = 1064, iterations = 1000, meanfun
## = median, : The following sampling units were removed because they were
## below sampling depth: 20RRNJune2017ITS, 25LRXJune2017ITS, 27RXXAugust2016ITS,
## 27RXXSeptember2017ITS, 32LXNDecember2017ITS, 46RXXMarch2017ITS,
## 47RRXDecember2017ITS, 48LXXSeptember2017ITS, 5LRNDecember2017ITS,
## 8RXNMarch2017ITS
bray.dist.ITS <- as.data.frame(as.matrix(median.avg.dist.ITS))
write.xlsx2(bray.dist.ITS, file = "bray.dist.ITS.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names = TRUE, append = FALSE)
NMDS.ITS <- metaMDS(bray.dist.ITS, distance = "bray", k = 3, noshare = FALSE, trymax = 400, trace = TRUE)
## Run 0 stress 0.1034416 
## Run 1 stress 0.102226 
## ... New best solution
## ... Procrustes: rmse 0.02181735  max resid 0.2444773 
## Run 2 stress 0.1034569 
## Run 3 stress 0.1052886 
## Run 4 stress 0.1022192 
## ... New best solution
## ... Procrustes: rmse 0.0106754  max resid 0.08277034 
## Run 5 stress 0.1066744 
## Run 6 stress 0.1022795 
## ... Procrustes: rmse 0.01290942  max resid 0.08527222 
## Run 7 stress 0.1021602 
## ... New best solution
## ... Procrustes: rmse 0.004342973  max resid 0.04197149 
## Run 8 stress 0.1023053 
## ... Procrustes: rmse 0.01315081  max resid 0.08762271 
## Run 9 stress 0.1018855 
## ... New best solution
## ... Procrustes: rmse 0.009392961  max resid 0.08142482 
## Run 10 stress 0.1049468 
## Run 11 stress 0.1024855 
## Run 12 stress 0.1022005 
## ... Procrustes: rmse 0.006606833  max resid 0.05844673 
## Run 13 stress 0.1028988 
## Run 14 stress 0.1034204 
## Run 15 stress 0.1020033 
## ... Procrustes: rmse 0.008730292  max resid 0.08205682 
## Run 16 stress 0.1018634 
## ... New best solution
## ... Procrustes: rmse 0.001273836  max resid 0.0101293 
## Run 17 stress 0.1035167 
## Run 18 stress 0.1036706 
## Run 19 stress 0.1021602 
## ... Procrustes: rmse 0.005574783  max resid 0.04095093 
## Run 20 stress 0.1028066 
## Run 21 stress 0.1047337 
## Run 22 stress 0.101726 
## ... New best solution
## ... Procrustes: rmse 0.007278107  max resid 0.08031383 
## Run 23 stress 0.1075901 
## Run 24 stress 0.1020344 
## ... Procrustes: rmse 0.008923795  max resid 0.0820247 
## Run 25 stress 0.103561 
## Run 26 stress 0.1087161 
## Run 27 stress 0.1094094 
## Run 28 stress 0.1050336 
## Run 29 stress 0.1036575 
## Run 30 stress 0.1049014 
## Run 31 stress 0.1033401 
## Run 32 stress 0.1046354 
## Run 33 stress 0.1026347 
## Run 34 stress 0.1034715 
## Run 35 stress 0.1069856 
## Run 36 stress 0.1045293 
## Run 37 stress 0.1032432 
## Run 38 stress 0.1036727 
## Run 39 stress 0.1017408 
## ... Procrustes: rmse 0.002186305  max resid 0.01568619 
## Run 40 stress 0.1022556 
## Run 41 stress 0.1033054 
## Run 42 stress 0.1047022 
## Run 43 stress 0.1025933 
## Run 44 stress 0.1057255 
## Run 45 stress 0.1016655 
## ... New best solution
## ... Procrustes: rmse 0.002711177  max resid 0.03583078 
## Run 46 stress 0.1070481 
## Run 47 stress 0.103016 
## Run 48 stress 0.1025487 
## Run 49 stress 0.1018629 
## ... Procrustes: rmse 0.007001865  max resid 0.07879939 
## Run 50 stress 0.1020274 
## ... Procrustes: rmse 0.008085198  max resid 0.08109908 
## Run 51 stress 0.1071755 
## Run 52 stress 0.1036474 
## Run 53 stress 0.1031518 
## Run 54 stress 0.1022378 
## Run 55 stress 0.1022273 
## Run 56 stress 0.1025939 
## Run 57 stress 0.1021349 
## ... Procrustes: rmse 0.008153416  max resid 0.07793045 
## Run 58 stress 0.1034201 
## Run 59 stress 0.1040976 
## Run 60 stress 0.1016764 
## ... Procrustes: rmse 0.001591751  max resid 0.017169 
## Run 61 stress 0.1036665 
## Run 62 stress 0.1108792 
## Run 63 stress 0.1030049 
## Run 64 stress 0.1026753 
## Run 65 stress 0.102592 
## Run 66 stress 0.1020972 
## ... Procrustes: rmse 0.006333918  max resid 0.06260569 
## Run 67 stress 0.1019168 
## ... Procrustes: rmse 0.007307465  max resid 0.08134799 
## Run 68 stress 0.1128113 
## Run 69 stress 0.1025593 
## Run 70 stress 0.1054687 
## Run 71 stress 0.109464 
## Run 72 stress 0.1016627 
## ... New best solution
## ... Procrustes: rmse 0.002280117  max resid 0.01837116 
## Run 73 stress 0.102266 
## Run 74 stress 0.1032318 
## Run 75 stress 0.1019791 
## ... Procrustes: rmse 0.008064308  max resid 0.08507382 
## Run 76 stress 0.102034 
## ... Procrustes: rmse 0.00815858  max resid 0.08194238 
## Run 77 stress 0.1030287 
## Run 78 stress 0.1045602 
## Run 79 stress 0.1037836 
## Run 80 stress 0.1052211 
## Run 81 stress 0.1027585 
## Run 82 stress 0.1029487 
## Run 83 stress 0.1025479 
## Run 84 stress 0.1025673 
## Run 85 stress 0.1017324 
## ... Procrustes: rmse 0.003368771  max resid 0.03510527 
## Run 86 stress 0.1022964 
## Run 87 stress 0.1036804 
## Run 88 stress 0.108209 
## Run 89 stress 0.1050465 
## Run 90 stress 0.102855 
## Run 91 stress 0.1020445 
## ... Procrustes: rmse 0.007367707  max resid 0.06864915 
## Run 92 stress 0.1024118 
## Run 93 stress 0.1027368 
## Run 94 stress 0.1055402 
## Run 95 stress 0.1043254 
## Run 96 stress 0.1052573 
## Run 97 stress 0.1023046 
## Run 98 stress 0.1057518 
## Run 99 stress 0.1044599 
## Run 100 stress 0.1049851 
## Run 101 stress 0.1038992 
## Run 102 stress 0.1025777 
## Run 103 stress 0.1086398 
## Run 104 stress 0.1052331 
## Run 105 stress 0.1037097 
## Run 106 stress 0.1017928 
## ... Procrustes: rmse 0.003752577  max resid 0.0350888 
## Run 107 stress 0.1053459 
## Run 108 stress 0.1030664 
## Run 109 stress 0.1026997 
## Run 110 stress 0.1036137 
## Run 111 stress 0.1019912 
## ... Procrustes: rmse 0.004516958  max resid 0.04281677 
## Run 112 stress 0.1023607 
## Run 113 stress 0.1053249 
## Run 114 stress 0.1018722 
## ... Procrustes: rmse 0.007185418  max resid 0.0827704 
## Run 115 stress 0.1049866 
## Run 116 stress 0.1027422 
## Run 117 stress 0.1016699 
## ... Procrustes: rmse 0.002573967  max resid 0.01607748 
## Run 118 stress 0.1032658 
## Run 119 stress 0.1038003 
## Run 120 stress 0.1037894 
## Run 121 stress 0.1028438 
## Run 122 stress 0.1022111 
## Run 123 stress 0.107964 
## Run 124 stress 0.1056293 
## Run 125 stress 0.1020346 
## ... Procrustes: rmse 0.007162788  max resid 0.06413547 
## Run 126 stress 0.1033928 
## Run 127 stress 0.1036479 
## Run 128 stress 0.1045988 
## Run 129 stress 0.1050291 
## Run 130 stress 0.1046954 
## Run 131 stress 0.1020148 
## ... Procrustes: rmse 0.005192132  max resid 0.04981483 
## Run 132 stress 0.1022318 
## Run 133 stress 0.1021045 
## ... Procrustes: rmse 0.00879418  max resid 0.08543031 
## Run 134 stress 0.103443 
## Run 135 stress 0.1033652 
## Run 136 stress 0.104009 
## Run 137 stress 0.1034414 
## Run 138 stress 0.1026358 
## Run 139 stress 0.1027934 
## Run 140 stress 0.103136 
## Run 141 stress 0.1045562 
## Run 142 stress 0.104973 
## Run 143 stress 0.1022534 
## Run 144 stress 0.1032806 
## Run 145 stress 0.10432 
## Run 146 stress 0.1025979 
## Run 147 stress 0.1025719 
## Run 148 stress 0.1147042 
## Run 149 stress 0.1027114 
## Run 150 stress 0.1052497 
## Run 151 stress 0.1026784 
## Run 152 stress 0.1043505 
## Run 153 stress 0.1025631 
## Run 154 stress 0.1035082 
## Run 155 stress 0.1047643 
## Run 156 stress 0.1025771 
## Run 157 stress 0.1020173 
## ... Procrustes: rmse 0.007927123  max resid 0.08239197 
## Run 158 stress 0.1026762 
## Run 159 stress 0.105241 
## Run 160 stress 0.1024549 
## Run 161 stress 0.1019162 
## ... Procrustes: rmse 0.005777645  max resid 0.06966694 
## Run 162 stress 0.1022609 
## Run 163 stress 0.1138861 
## Run 164 stress 0.1041216 
## Run 165 stress 0.1017235 
## ... Procrustes: rmse 0.003172806  max resid 0.03560678 
## Run 166 stress 0.1039789 
## Run 167 stress 0.1025309 
## Run 168 stress 0.1019256 
## ... Procrustes: rmse 0.005060735  max resid 0.05253533 
## Run 169 stress 0.1033446 
## Run 170 stress 0.1139195 
## Run 171 stress 0.1030149 
## Run 172 stress 0.1018958 
## ... Procrustes: rmse 0.005181505  max resid 0.05490719 
## Run 173 stress 0.1051713 
## Run 174 stress 0.1035176 
## Run 175 stress 0.103605 
## Run 176 stress 0.10508 
## Run 177 stress 0.1030392 
## Run 178 stress 0.1052827 
## Run 179 stress 0.1027528 
## Run 180 stress 0.1023023 
## Run 181 stress 0.1018591 
## ... Procrustes: rmse 0.004564508  max resid 0.05245591 
## Run 182 stress 0.1017523 
## ... Procrustes: rmse 0.004188954  max resid 0.03464671 
## Run 183 stress 0.1072388 
## Run 184 stress 0.1033439 
## Run 185 stress 0.1025586 
## Run 186 stress 0.1029229 
## Run 187 stress 0.1025425 
## Run 188 stress 0.1020394 
## ... Procrustes: rmse 0.007627395  max resid 0.06618612 
## Run 189 stress 0.1018584 
## ... Procrustes: rmse 0.006902274  max resid 0.08236374 
## Run 190 stress 0.1030269 
## Run 191 stress 0.1029932 
## Run 192 stress 0.10464 
## Run 193 stress 0.1038516 
## Run 194 stress 0.1018513 
## ... Procrustes: rmse 0.006737391  max resid 0.08211791 
## Run 195 stress 0.1025763 
## Run 196 stress 0.1155397 
## Run 197 stress 0.1053756 
## Run 198 stress 0.1025485 
## Run 199 stress 0.1033358 
## Run 200 stress 0.105229 
## Run 201 stress 0.1038723 
## Run 202 stress 0.1019398 
## ... Procrustes: rmse 0.005902901  max resid 0.05628087 
## Run 203 stress 0.1034608 
## Run 204 stress 0.1028819 
## Run 205 stress 0.1053369 
## Run 206 stress 0.1035027 
## Run 207 stress 0.1051262 
## Run 208 stress 0.1075342 
## Run 209 stress 0.1026542 
## Run 210 stress 0.1074899 
## Run 211 stress 0.1045322 
## Run 212 stress 0.1050699 
## Run 213 stress 0.1032503 
## Run 214 stress 0.1026424 
## Run 215 stress 0.1037207 
## Run 216 stress 0.1070578 
## Run 217 stress 0.1023244 
## Run 218 stress 0.1023371 
## Run 219 stress 0.1081403 
## Run 220 stress 0.1025895 
## Run 221 stress 0.1038863 
## Run 222 stress 0.1035582 
## Run 223 stress 0.1061032 
## Run 224 stress 0.1023903 
## Run 225 stress 0.1027721 
## Run 226 stress 0.1025444 
## Run 227 stress 0.1051926 
## Run 228 stress 0.1032486 
## Run 229 stress 0.1025692 
## Run 230 stress 0.1033395 
## Run 231 stress 0.1027101 
## Run 232 stress 0.1036359 
## Run 233 stress 0.1033371 
## Run 234 stress 0.1017805 
## ... Procrustes: rmse 0.004208538  max resid 0.04196572 
## Run 235 stress 0.1020692 
## ... Procrustes: rmse 0.007853111  max resid 0.07078279 
## Run 236 stress 0.1027522 
## Run 237 stress 0.1022305 
## Run 238 stress 0.1035704 
## Run 239 stress 0.1073083 
## Run 240 stress 0.1034984 
## Run 241 stress 0.1052264 
## Run 242 stress 0.1019462 
## ... Procrustes: rmse 0.007695539  max resid 0.08465692 
## Run 243 stress 0.1026462 
## Run 244 stress 0.1034169 
## Run 245 stress 0.1024427 
## Run 246 stress 0.1057999 
## Run 247 stress 0.103368 
## Run 248 stress 0.1043746 
## Run 249 stress 0.1025407 
## Run 250 stress 0.1018606 
## ... Procrustes: rmse 0.004734909  max resid 0.05428239 
## Run 251 stress 0.1046786 
## Run 252 stress 0.1024555 
## Run 253 stress 0.1055564 
## Run 254 stress 0.1019498 
## ... Procrustes: rmse 0.005515581  max resid 0.05457991 
## Run 255 stress 0.101737 
## ... Procrustes: rmse 0.00339964  max resid 0.03563156 
## Run 256 stress 0.1020452 
## ... Procrustes: rmse 0.00844098  max resid 0.08612605 
## Run 257 stress 0.1045068 
## Run 258 stress 0.1017399 
## ... Procrustes: rmse 0.004131214  max resid 0.02551933 
## Run 259 stress 0.1022028 
## Run 260 stress 0.101889 
## ... Procrustes: rmse 0.00666658  max resid 0.07969948 
## Run 261 stress 0.1043615 
## Run 262 stress 0.1024206 
## Run 263 stress 0.1052071 
## Run 264 stress 0.1045263 
## Run 265 stress 0.1020916 
## ... Procrustes: rmse 0.008594322  max resid 0.07611796 
## Run 266 stress 0.1049277 
## Run 267 stress 0.1049517 
## Run 268 stress 0.1054611 
## Run 269 stress 0.1032181 
## Run 270 stress 0.1047857 
## Run 271 stress 0.1045819 
## Run 272 stress 0.103215 
## Run 273 stress 0.1034907 
## Run 274 stress 0.105452 
## Run 275 stress 0.1045778 
## Run 276 stress 0.103836 
## Run 277 stress 0.1054566 
## Run 278 stress 0.1051983 
## Run 279 stress 0.102321 
## Run 280 stress 0.1035903 
## Run 281 stress 0.1033324 
## Run 282 stress 0.1016708 
## ... Procrustes: rmse 0.002194432  max resid 0.02447812 
## Run 283 stress 0.1045206 
## Run 284 stress 0.1064861 
## Run 285 stress 0.1037251 
## Run 286 stress 0.1074868 
## Run 287 stress 0.1027372 
## Run 288 stress 0.1022857 
## Run 289 stress 0.1050393 
## Run 290 stress 0.1090944 
## Run 291 stress 0.1037347 
## Run 292 stress 0.1031956 
## Run 293 stress 0.1049937 
## Run 294 stress 0.103446 
## Run 295 stress 0.1016554 
## ... New best solution
## ... Procrustes: rmse 0.001915193  max resid 0.01460278 
## Run 296 stress 0.1026395 
## Run 297 stress 0.1032238 
## Run 298 stress 0.108303 
## Run 299 stress 0.105138 
## Run 300 stress 0.1056164 
## Run 301 stress 0.1045625 
## Run 302 stress 0.1027069 
## Run 303 stress 0.103958 
## Run 304 stress 0.1025547 
## Run 305 stress 0.1037276 
## Run 306 stress 0.1017485 
## ... Procrustes: rmse 0.003497165  max resid 0.02062046 
## Run 307 stress 0.1019178 
## ... Procrustes: rmse 0.004790139  max resid 0.05411592 
## Run 308 stress 0.1020442 
## ... Procrustes: rmse 0.005590636  max resid 0.05534946 
## Run 309 stress 0.1060374 
## Run 310 stress 0.1026752 
## Run 311 stress 0.1018531 
## ... Procrustes: rmse 0.007001605  max resid 0.08187253 
## Run 312 stress 0.106804 
## Run 313 stress 0.1044966 
## Run 314 stress 0.1020695 
## ... Procrustes: rmse 0.005723969  max resid 0.05691739 
## Run 315 stress 0.1036345 
## Run 316 stress 0.103325 
## Run 317 stress 0.1019269 
## ... Procrustes: rmse 0.004949187  max resid 0.05580597 
## Run 318 stress 0.1022977 
## Run 319 stress 0.1026207 
## Run 320 stress 0.1025068 
## Run 321 stress 0.1034426 
## Run 322 stress 0.1037351 
## Run 323 stress 0.1060372 
## Run 324 stress 0.1020177 
## ... Procrustes: rmse 0.007898818  max resid 0.08071515 
## Run 325 stress 0.1035978 
## Run 326 stress 0.1026315 
## Run 327 stress 0.1022312 
## Run 328 stress 0.1026271 
## Run 329 stress 0.1028372 
## Run 330 stress 0.1029809 
## Run 331 stress 0.1071241 
## Run 332 stress 0.1078917 
## Run 333 stress 0.1025533 
## Run 334 stress 0.104177 
## Run 335 stress 0.1044862 
## Run 336 stress 0.1050522 
## Run 337 stress 0.1016668 
## ... Procrustes: rmse 0.002204927  max resid 0.02284503 
## Run 338 stress 0.1035263 
## Run 339 stress 0.1050347 
## Run 340 stress 0.105228 
## Run 341 stress 0.1018703 
## ... Procrustes: rmse 0.004501693  max resid 0.05519265 
## Run 342 stress 0.1054364 
## Run 343 stress 0.102541 
## Run 344 stress 0.1043361 
## Run 345 stress 0.1029396 
## Run 346 stress 0.1049648 
## Run 347 stress 0.1027477 
## Run 348 stress 0.1034875 
## Run 349 stress 0.1051458 
## Run 350 stress 0.103004 
## Run 351 stress 0.1040485 
## Run 352 stress 0.1035455 
## Run 353 stress 0.1022223 
## Run 354 stress 0.1017381 
## ... Procrustes: rmse 0.003034046  max resid 0.03617134 
## Run 355 stress 0.1025491 
## Run 356 stress 0.1108277 
## Run 357 stress 0.1021067 
## ... Procrustes: rmse 0.006369138  max resid 0.05852839 
## Run 358 stress 0.1064937 
## Run 359 stress 0.1088592 
## Run 360 stress 0.1051294 
## Run 361 stress 0.1018992 
## ... Procrustes: rmse 0.004295914  max resid 0.05028761 
## Run 362 stress 0.1022957 
## Run 363 stress 0.1023214 
## Run 364 stress 0.1062385 
## Run 365 stress 0.1033258 
## Run 366 stress 0.1027441 
## Run 367 stress 0.1052934 
## Run 368 stress 0.101762 
## ... Procrustes: rmse 0.00351508  max resid 0.03640896 
## Run 369 stress 0.1021703 
## Run 370 stress 0.1030053 
## Run 371 stress 0.103268 
## Run 372 stress 0.1050429 
## Run 373 stress 0.1035997 
## Run 374 stress 0.103363 
## Run 375 stress 0.1020397 
## ... Procrustes: rmse 0.007650374  max resid 0.08010203 
## Run 376 stress 0.1057458 
## Run 377 stress 0.1031281 
## Run 378 stress 0.1018803 
## ... Procrustes: rmse 0.004253365  max resid 0.05252094 
## Run 379 stress 0.1026492 
## Run 380 stress 0.1049269 
## Run 381 stress 0.102607 
## Run 382 stress 0.1034051 
## Run 383 stress 0.1022014 
## Run 384 stress 0.104962 
## Run 385 stress 0.1059308 
## Run 386 stress 0.103671 
## Run 387 stress 0.1021671 
## Run 388 stress 0.1025342 
## Run 389 stress 0.1036991 
## Run 390 stress 0.1026817 
## Run 391 stress 0.1037019 
## Run 392 stress 0.101951 
## ... Procrustes: rmse 0.005781329  max resid 0.05641397 
## Run 393 stress 0.102302 
## Run 394 stress 0.1099411 
## Run 395 stress 0.1031918 
## Run 396 stress 0.1034833 
## Run 397 stress 0.1027881 
## Run 398 stress 0.1054689 
## Run 399 stress 0.1028868 
## Run 400 stress 0.1059391 
## *** No convergence -- monoMDS stopping criteria:
##    291: no. of iterations >= maxit
##    109: stress ratio > sratmax
stressplot(NMDS.ITS)

coordinates.ITS <- data.frame(NMDS.ITS$points[,1:3])
nmds.metadata.ITS <- merge(coordinates.ITS, metadata.ITS.1, by.x = "row.names", by.y = "row.names")
nmds.metadata.ITS <- data.frame(nmds.metadata.ITS, col.names    =   TRUE,   row.names   =   TRUE)

fungi.nmds.plot <-  ggplot(data = nmds.metadata.ITS) +
                       aes(x = MDS1, y = MDS3, color = as.factor(Precipitation)) +
                       geom_point(aes(shape = as.factor(Ecosystem), size = as.factor(Ecosystem))) +
                       scale_color_manual(values = c('#408080', '#FF8040')) +
                       scale_shape_manual(values = c(17,16)) +
                       scale_size_manual(values = c(5,5)) +
                       geom_text(label = as.factor(nmds.metadata.ITS$eco.prec), size = 1) +
                       theme_minimal() + 
                       theme(legend.position = "none") +
                       labs(col = "Precipitation", shape = "Ecosystem") + 
                       theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                       panel.background = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=1), axis.text.x=element_blank(), 
                       axis.ticks.x=element_blank(),axis.text.y=element_blank(),  
                       axis.ticks.y=element_blank())

plot(fungi.nmds.plot)

PERMANOVA.ITS <- adonis(bray.dist.ITS ~ Block + Ecosystem*Precipitation + Ecosystem*Nitrogen + Ecosystem*Collection_date + Collection_date*Precipitation + Collection_date*Nitrogen + Precipitation*Nitrogen, data = nmds.metadata.ITS, strata = nmds.metadata.ITS$Block, permutations = 999)

PERMANOVA.ITS
## 
## Call:
## adonis(formula = bray.dist.ITS ~ Block + Ecosystem * Precipitation +      Ecosystem * Nitrogen + Ecosystem * Collection_date + Collection_date *      Precipitation + Collection_date * Nitrogen + Precipitation *      Nitrogen, data = nmds.metadata.ITS, permutations = 999, strata = nmds.metadata.ITS$Block) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                                Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## Block                           1     1.657 1.65720  7.2005 0.03203  0.001 ***
## Ecosystem                       1     1.074 1.07426  4.6676 0.02076  0.001 ***
## Precipitation                   1     0.615 0.61529  2.6734 0.01189  0.002 ** 
## Nitrogen                        1     0.213 0.21280  0.9246 0.00411  0.495    
## Collection_date                 6     2.435 0.40575  1.7630 0.04705  0.001 ***
## Ecosystem:Precipitation         1     0.349 0.34851  1.5143 0.00674  0.070 .  
## Ecosystem:Nitrogen              1     0.209 0.20875  0.9070 0.00403  0.502    
## Ecosystem:Collection_date       6     2.225 0.37078  1.6110 0.04300  0.001 ***
## Precipitation:Collection_date   6     1.848 0.30796  1.3381 0.03571  0.022 *  
## Nitrogen:Collection_date        6     1.404 0.23400  1.0167 0.02714  0.399    
## Precipitation:Nitrogen          1     0.355 0.35459  1.5407 0.00685  0.075 .  
## Residuals                     171    39.356 0.23015         0.76067           
## Total                         202    51.738                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###Plant###
median.avg.dist.plant <- vegdist(otu.plant.2015.1, method = "bray")
bray.dist.plant <- as.data.frame(as.matrix(median.avg.dist.plant))
write.xlsx2(bray.dist.plant, file = "bray.dist.plant.xlsx", sheetName = "Sheet1", col.names = TRUE, row.names = TRUE, append = FALSE)
NMDS.plant <- metaMDS(bray.dist.plant, distance = "bray", k = 3, noshare = FALSE, trymax = 400, trace = TRUE)
## Run 0 stress 0.1019617 
## Run 1 stress 0.1022528 
## ... Procrustes: rmse 0.03165542  max resid 0.08401119 
## Run 2 stress 0.1019806 
## ... Procrustes: rmse 0.01193097  max resid 0.03951897 
## Run 3 stress 0.1003177 
## ... New best solution
## ... Procrustes: rmse 0.03954919  max resid 0.141736 
## Run 4 stress 0.1047205 
## Run 5 stress 0.1045815 
## Run 6 stress 0.1084601 
## Run 7 stress 0.100317 
## ... New best solution
## ... Procrustes: rmse 0.0005105433  max resid 0.001521973 
## ... Similar to previous best
## Run 8 stress 0.1045886 
## Run 9 stress 0.1019759 
## Run 10 stress 0.1003173 
## ... Procrustes: rmse 0.0003661874  max resid 0.001092883 
## ... Similar to previous best
## Run 11 stress 0.1084616 
## Run 12 stress 0.1003181 
## ... Procrustes: rmse 0.0004177886  max resid 0.001676165 
## ... Similar to previous best
## Run 13 stress 0.1052178 
## Run 14 stress 0.1019863 
## Run 15 stress 0.1019478 
## Run 16 stress 0.1093034 
## Run 17 stress 0.1084781 
## Run 18 stress 0.1019414 
## Run 19 stress 0.101944 
## Run 20 stress 0.1003173 
## ... Procrustes: rmse 0.0003675437  max resid 0.001147174 
## ... Similar to previous best
## *** Solution reached
stressplot(NMDS.plant)

coordinates.plant <- data.frame(NMDS.plant$points[,1:3])
nmds.metadata.plant <- merge(coordinates.plant, metadata.plant.2015.1, by.x = "row.names", by.y = "row.names")
nmds.metadata.plant <- data.frame(nmds.metadata.plant, col.names    =   TRUE,   row.names   =   TRUE)

plant.nmds.plot <- ggplot(data = nmds.metadata.plant) +
                   aes(x = MDS1, y = MDS2, color = as.factor(precipitation)) + 
                   geom_point(aes(shape = as.factor(ecosystem), size = as.factor(ecosystem))) +
                   scale_color_manual(values = c('#408080', '#FF8040')) +
                   scale_shape_manual(values = c(17,16)) +
                   scale_size_manual(values = c(5,5)) +
                   geom_text(label = as.factor(nmds.metadata.plant$eco.prec), size = 1) +
                   theme_minimal() + 
                   theme(legend.position = "right") +
                   labs(col = "Precipitation", shape = "Ecosystem") + 
                   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                   panel.background = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=1), axis.text.x=element_blank(), 
                   axis.ticks.x=element_blank(),axis.text.y=element_blank(), 
                   axis.ticks.y=element_blank())

plot(plant.nmds.plot)

bray.dist.plant <- bray.dist.plant[(names(bray.dist.plant) %in% row.names(metadata.plant.2015.1)),]
bray.dist.plant <- bray.dist.plant[,(names(bray.dist.plant) %in% row.names(bray.dist.plant))]
matching.metadata.bray.plant <- metadata.plant.2015.1[c(match(sort(row.names(bray.dist.plant)), sort(row.names(metadata.plant.2015.1)))),]

PERMANOVA.plant <- adonis(bray.dist.plant ~ block + ecosystem*precipitation  + precipitation*nitrogen + ecosystem*nitrogen, data = matching.metadata.bray.plant, strata = matching.metadata.bray.plant$block,  permutations = 999)

PERMANOVA.plant
## 
## Call:
## adonis(formula = bray.dist.plant ~ block + ecosystem * precipitation +      precipitation * nitrogen + ecosystem * nitrogen, data = matching.metadata.bray.plant,      permutations = 999, strata = matching.metadata.bray.plant$block) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                         Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## block                    1    2.2688 2.26885 18.6187 0.26810  0.001 ***
## ecosystem                1    0.9055 0.90554  7.4311 0.10700  0.001 ***
## precipitation            1    1.4921 1.49207 12.2443 0.17631  0.001 ***
## nitrogen                 1    0.0778 0.07777  0.6382 0.00919  0.659    
## ecosystem:precipitation  1    0.6702 0.67015  5.4994 0.07919  0.006 ** 
## precipitation:nitrogen   1    0.0565 0.05651  0.4638 0.00668  0.795    
## ecosystem:nitrogen       1    0.0672 0.06716  0.5511 0.00794  0.686    
## Residuals               24    2.9246 0.12186         0.34559           
## Total                   31    8.4626                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#stacked bar plot of variance by factors
variance <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/variance_fig.xlsx")
variance <- as.data.frame(variance, col.names   =   TRUE,   row.names   =   TRUE)
variance[1:3 , 1:3]
##      factor community  percent
## 1 Residuals  Bacteria 59.86696
## 2 Residuals     Fungi 81.81324
## 3 Residuals     Plant 22.38277
c11 <- c("#cacccc",
         "#4f4447",
         "#697456",
         "#91ab98",
         "#504254",
         "#408080",
         "#54bfab",
         "#306773",
         "#809a97",
         "#b48a6f",
         "#9d7e72")

variance.plot <- ggplot(variance, aes(x = community, y = percent))+
                 geom_col(aes(fill = factor(factor, levels = c("Residuals", "Block", "Ecosystem", "Ecosystem_collection date","Collection_Date", "Precipitation", "Ecosystem_precipitation", 
                                                               "Precipitation_collection date", "Precipitation_nitrogen", "Nitrogen", "Ecosystem_nitrogen"))), width = 0.7) +
                 theme_minimal() +  
                 theme(axis.line.x=element_line(colour = "black"), axis.line.y=element_line(colour = "black")) + 
                 scale_fill_manual(values = c11) +
                 labs(fill = "Factor (s)") +
                 #ggtitle("Variance Explained By Factor") + #Adds tittle and subtittle. Can modify p and r-squared values + title.
                 guides(size = "none") # makes the sizes guide dissapear
  
plot(variance.plot)

combined.nmds.variance.plots <- ggarrange(bacteria.nmds.plot, fungi.nmds.plot, plant.nmds.plot, variance.plot + font("x.text", size = 10), labels = c("A.", "B.", "C.", "D."), ncol = 2, nrow = 2, widths = c(5,5,5,4))
plot(combined.nmds.variance.plots)

ambient.water.input <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/ambient_water_input_figure.xlsx")
ambient.water.input <- as.data.frame(ambient.water.input, col.names =   TRUE,   row.names   =   TRUE)
ambient.water.input[1:2 , 1:3]
##   month_year       date ambient_water_input_mm
## 1   oct_2014 2014-10-01                      0
## 2   oct_2014 2014-10-02                      0
sapply(ambient.water.input, class)
## $month_year
## [1] "character"
## 
## $date
## [1] "POSIXct" "POSIXt" 
## 
## $ambient_water_input_mm
## [1] "numeric"
ambient.water.input.1 <- ambient.water.input %>% mutate_if(is.POSIXt, as.Date)
ambient.water.input.2 <- ambient.water.input.1 %>% mutate_if(is.POSIXct, as.Date)
sapply(ambient.water.input.2, class)
##             month_year                   date ambient_water_input_mm 
##            "character"                 "Date"              "numeric"
reduced.water.input <- read_excel("~/Google_drive/my_projects/lrgce/final_analysis_finks_etal_2020/reduced_water_input_figure.xlsx")
reduced.water.input <- as.data.frame(reduced.water.input, col.names =   TRUE,   row.names   =   TRUE)
reduced.water.input[1:2 , 1:3]
##   month_year       date reduced_water_input_mm
## 1   oct_2014 2014-10-01                      0
## 2   oct_2014 2014-10-02                      0
sapply(reduced.water.input, class)
## $month_year
## [1] "character"
## 
## $date
## [1] "POSIXct" "POSIXt" 
## 
## $reduced_water_input_mm
## [1] "numeric"
reduced.water.input.1 <- reduced.water.input %>% mutate_if(is.POSIXt, as.Date)
reduced.water.input.2 <- reduced.water.input.1 %>% mutate_if(is.POSIXct, as.Date)
sapply(reduced.water.input.2, class)
##             month_year                   date reduced_water_input_mm 
##            "character"                 "Date"              "numeric"
pd = position_dodge(0.5)

water.fig <- ggplot() + geom_point(data = ambient.water.input.2, aes(x = date, y = ambient_water_input_mm), size=2, pch = 19, color = "#408080") +
             geom_point(data = reduced.water.input.2, aes(x = date, y = reduced_water_input_mm), size=2, pch = 19, position = pd, color = "#a9805d") +
             xlab('Dates') +
             ylab('Precipitation (mm)') +
             geom_line(data = ambient.water.input.2, aes(x = date, y = ambient_water_input_mm), color = "#408080") +
             geom_line(data = reduced.water.input.2, aes(x = date, y = reduced_water_input_mm), position = pd, color = "#a9805d") +
             scale_y_continuous(breaks = seq(0,70,by=5), expand = c(0, 0)) +
             scale_x_date(breaks = date_breaks("month"), labels = date_format("%b"), expand = c(0, 0)) +
             theme(panel.grid.major = element_line(colour = "black", linetype = "dotted"), panel.grid.minor = element_blank(), panel.grid.major.x= element_blank(),
             panel.background = element_blank(), panel.border = element_rect(colour = "black", fill=NA, size=0.5)) +
             scale_linetype_manual(values = c(1,1))
  
plot(water.fig)
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).